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REPORT  SUMMARY 


This  final  report  covers  the  nine  month  period  from 
July  1976  through  the  termination  of  the  contract  at  the  end 
of  March  1977.  As  such,  this  report  describes  the 
accomplishments  of  this  period  in  detail  and  thus  includes 
the  semi-annual  report  of  the  period  July  1976  - December 
1976  . 


Section  I discusses  the  theoretical  underpinnings  of  a 
new  mathematical  model  describing  some  of  the  perceptual 
characteristics  of  human  sensory  information  processing. 
While  much  of  the  mathematics  of  the  proposed  model  is 
explained  here,  experimental  testing  and  verification  has 
yet  to  take  place.  This  work  will  continue  to  be  pursued 
under  alternate  funding. 

Section  II  discusses  in  detail  the  background 
considerations  necessary  to  understand  the  problems  of 
removing  atmospheric  turbulence  blur  from  images.  This  area 
of  research  is  currently  in  transition  from  the  stage  of 
artificial  computer  models,  to  the  stage  of  working  with  the 
simplest  real  data,  i.e.  case  of  a point  light  source,  a 
star,  viewed  through  the  earth's  atmosphere.  After  the 
termination  of  this  contract  this  work  will  be  funded  by  the 
Air  Force  Office  of  Scientific  Research  under  grant  number 
77-3212. 
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Section  III  reports  the  background  and  progress  to  date 
on  building  the  tools  necessary  for  an  image  understanding 
system.  These  tools  as  described  were  incorporated  into  a 
manual  image  understanding  system  designed  to  do  analysis  by 
synthesis.  The  system  handles  two-dimensional  single 
objects  on  non-textured  backgrounds,  but  only  at  the  level 
of  a trial  testbed.  Necessary  next  steps  in  developing  a 
working  system  are  outlined. 
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SECTION  I 

GROUP  REPRESENTATIONS  AND  THE  MODELING 

OF  VISUAL  PERCEPTION 


James  T.  Kajiya 


The  student  of  human  visual  perception  is  often 
overwhelmed  by  the  vast  amount  of  data  that  has  been 
ac cumuli  ted  from  experiments  performed  within  the  last 
century  or  so.  It  is  often  difficult  to  understand  why  a 
certain  experiment  has  been  performed.  Results  from  similar 
experiments  sometimes  seem  to  conflict.  Further  confusion 
results  when  the  student  encounters  raging  controversies, 
the  resolution  of  which  would  seem  to  minimally  advance  our 
knowledge  of  how  we  see.  The  reason  for  all  this  trouble 
stems  from  the  fact  that  a suitable  superstructure  providing 
organization  and  support  of  this  accumulation  of  data  does 
not  exist,  i.e.  an  adequate  theory  of  perception  is  not  in 
hand  . 


Theories,  by  their  nature,  are  primarily  qualitative 
descriptions  of  what  is  going  on.  Of  the  qialitative 
observations  that  can  be  made  about  the  human  visual  system, 
some  of  the  most  fruitful  are  those  of  the  perceptual 
constancies.  In  short,  psychophysi cists  have  puzzled  over 
the  eye's  marvelous  capacity  to  respond  gracefully  over  wide 
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transformations  on  the  images  presented  to  the  retina.  The 
perception  of  a triangle  is  relatively  unaffected  by 
variations  in  size,  orientation,  and  perspective.  This 
particular  phenomena  takes  as  its  name  "shape  constancy". 
These  constancies  or  invariances  provide  enough  qualitative 
information  to  draw  upon  the  powerful  theory  of  group 
representations. 

This  paper  will  report  on  some  first  steps  taken  toward 
the  application  of  the  theory  of  group  representations  to 
the  modeling  of  human  perception.  In  particular,  Section  2 
will  briefly  discuss  two  groups  that  express  certain 
observed  invariances  within  the  visual  system.  Section  3 
will  perform  the  relatively  straightforward  calculations  of 
the  unitary  equivalence  classes  of  irreducible 
representations  for  one  of  these  groups.  Section  4 will 
attempt  to  carry  out  a brute  force  harmonic  analysis.  It 
will  turn  out  that  this  brute  force  approach  will  fail  in 
tnat  it  will  not  shed  any  light  upon  the  workings  of  the 
visual  system.  This  may  explain  in  part  why  Abstract 
Harmonic  Analysis  has  not  seen  wide  use  in  the  study  of 
vision  to  date.  Section  5 (the  principal  contribution  of 
the  paper)  will  examine  why  the  brute  force  approach  fails. 
Also  an  operator  imbedding  the  image  space  into  a larger 
3pace  is  presented  for  the  first  time.  This  new  operator 
will  be  seen  to  possess  certain  attractive  properties  that 
may  turn  out  to  make  it  quite  useful  as  a gadget  for  image 
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processing  in  general.  Finally,  Section  6 will  quickly 


outline  the  rest  of  the  program  needed  to  obtain 


perceptually  important  transform  that  will  mimic  the 


operation  of  the  visual  system. 


2.  THE  PERCEPTUAL  GROUPS 


It  has  been  known  for  quite  some  time  that  the  groups 


describing  the  various  perceptual  constancies  are  Lie  Groups 


tl)  2],  In  particular,  Pitts  and  McCulloch  in  [1]  conducted 


a kind  of  harmonic  analysis  at  the  "DC"  frequency.  The  main 


group  that  should  be  kept  in  mind  is  the  group  G of  affine 


transformations  on  8.  Our  group  G expresses  the  form 


constancy  in  the  literature.  Let  (x,y,w)  be  homogeneous 


co-ordinates  of  the  plane.  The  one  parameter  subgroups 


generating  G are: 


translations 


1 0 a 


0 1 b 


0 0 ? 


rotat ions 


cose  -sine  0 


sine  cose 


dilatations 


I . 


reflection 


Thir-  group  is  also  a serai-direct  product  of  two  abelian 
groups.  In  what  follows  we  treat  only  the  first  group. 


Notice  that  if  the  identification  of  8^  with  the  complex 
plane  is  made  then  G is  the  complex  "otx+B"  group  where 
3 = (a,b) 


Since  the  first  derived  subalgebra  >(^  ] of  the  Lie 

Algebra  (JJ  of  G is  the  set  of  translations  which  is  abelian, 
this  group  is  solvable.  It  is  also  apparent  from  the 


otx  + B 


expression  that  G is  the  semi-direct  product  of  the 


two  abelian  groups  (C,  + ) and  (I\{0)  ,•). 


There  is  some  reason  to  suspect  that  rotations  do  not 
correspond  to  a perceptual  invariance  [3,  4],  Furthermore 
bilateral  symmetry  seems  to  evoke  a strong  mechanism  in  the 
visual  system.  Indeed,  study  of  this  kind  of  symmetry  dates 
back  to  Mach.  Thus  another  important  group  we  should 
consider  is  generated  by  translation,  dilatation,  and 


m’r&mnrmim*  tm 
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Since  these  groups  have  such  a nice  structure,  it  is 
quite  straightforward  to  apply  the  Mackey  Induction  Theory 
to  obtain  a complete  characterization  of  their  unitary 

a 

equivalence  classes  of  irreducible  representations  G. 

3.  CALCULATION  OF  THE  DUAL  6 


Let  us  recall  some  basic  definitions  and  theorems  [5, 
6].  Definition.  Let  Y be  a representation  of  a subgroup  H 


of  G on  a hilbert  space^y  Consider  the  space^  of  all 


functions  from  G to»^ y satisfying 


i) 
i i ) 
i i i ) 


x (f  (x)  ,5  ) is  Borel 


f(hx)  = y(h)f(x) , heH 
f I |f(x)| |2  dx  < 


JG/H 


almost  al  1 


oo 


xeG 


where  we  make  the  usual  identification.  Set 


7r(g)f(x)  = f(xg)  (p)1/2  f 


and  call  it  the  representation  of  G induced  from  H by  y,  We 

Q 

write  tt  = I n d ^ ( y ) . 

Theorem . Let  G be  a regular  semi-direct  product  of  H and  N, 

A 

with  N abelian  and  normal.  For  each  orbit  H *Y,  YeN  , set 

A 

H ={heH:h*Y=Y}  Get  u ( x ,y ) = VxWyeG  ' Then  the 

A 

projection-valued  measure  defined  by  W in  N is  concentrated 


w 
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U 

in  a single  orbit  H*v  and  V = IndM  (v)  for  some  veH  . 

H y 

Y 

Every  pair  consisting  of  an  orbit  H*y  and  irreducible  V 
arise  in  this  way,  and  two  of  these  are  equivalent  if  the 
orbits  are  identical  and  the  representations  \>  are 
equivalent. 

The  above  definition  and  theory  were  taken  from  the 
excellent  [ 7 ] . 

Since  both  perceptual  groups  have  available  the  abelian 
2 

normal  subgroup  1 , the  characters  in  question  y are  the 

2 

familiar  fourier  kernels  which  are  isomorphic  to  tt  . The 
orbits  form  two  disjoint  subsets  of  the  plane:  The  point 
(0,0)  and  the  punctured  plane.  Their  respective  isotropy 
subgroups  are  H itself  and  the  identity. 

4.  A NAIVE  ATTEMPT 


Since  we  are  interested  in  images,  the  hilbert  space 
that  our  representation  acts  upon  will  be  sl  o(l^,  p ) 
The  action  of  G on  f)  will  be 


ug(f(x))  = a(g)f(gx) 

The  multiplier  a(g)  is  the  square  root  of  the  Radon-Nikodym 
derivative  of  t,  he  measure 


P g ( E ) = y ( g E ) with  respect  to 
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Lesbegue  measure  y.  It  is  necessary  because,  as  a moments 
reflection  will  show,  the  best  we  can  do  is  quasi-invariance 
of  our  measure.  The  a of  course  makes  U a unitary 
representation  . 

Our  aim  is  to  decompose  this  representation  into 
irreducible  subrepresentations  and  then  reap  the  benefits  by 
identifying  the  projection-valued  measures.  Unfortunately 
this  representation  is  itself  irreducible  so  our  hopes  of 
arriving  at  a decomposition  have  evaporated.  The 

irreducibility  of  this  representation  indicates  that  any 
ad-hoc  scheme  that  attempts  to  "crank  in"  rotations  and 
dilatations  will  meet  with  only  partial  success  at  best. 
The  main  problem  with  our  space  of  images  is  that  there  is 
not  enough  "elbow  room"  to  be  able  to  proceed  with  our 
decomposition.  To  see  this  we  present  a "heuristic" 
argument.  Let  f(x)  be  the  characteristic  function  of,  say, 
the  unit  square.  We  can  approximate  any  characteristic 
function  of  a measurable  subset  of  the  plane  simply  by 
following  the  classical  Lesbegue  construction.  The  linear 
closure  of  these  functions  form  the  simple  functions.  Since 
the  simple  functions  are  dense  in  we  see  that  any  picture 
is  related  to  any  other. 


I | 
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5.  THE  MANDALA  TRANSFORM 

Obviously  something  must  be  wrong  with  our  model  of  the 
human  visual  system.  The  eye  can  successfully  decompose 
images  into  separate  parts.  The  discrepancy  can  be 
accounted  for  by  the  notion  of  texture . Key  to  the  above 
heuristic  argument  was  the  freedom  to  shrink  the  unit  square 
to  arbitrarily  small  proportions.  If  this  is  done  with  an 
image  then  at  some  point  the  image  will  change  its  perceived 
character.  Objects  in  the  image  will  become  so  sm~ll  that 
they  will  become  indiscernible  and  turn  into  texture.  A 
beautiful  example  of  this  effect  are  the  images  by  I irmon 
and  Knowlton  [8].  In  other  words  the  dilatation  operation 
is  not  a true  symmetry  operation.  How  then  can  we  have  form 
constancy? 

Further  clues  to  the  operation  of  the  visual  system  can 
be  culled  from  the  landmark  experiments  of  Hube]  and  Wiesel 
[9]  . In  these  experiments  the  single  neurons  in  the  visual 
cortex  monitored  by  microelectrodes  were  found  to  respond 
only  to  small  portions  of  the  visual  field.  An  examination 
of  several  closely  spaced  neurons  indicated  slightly 
displaced  but  overlapping  "windows"  into  the  visual  field. 
Furthermore,  many  neurons  were  responsive  only  to  lines 
placed  in  a certain  orientation  with  respect  to  the  window. 
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Hubei  and  Wiesel  concluded  that  they  had  found  mechanisms 
closely  resembling  "line  detectors"  within  the  visual 
cortex.  Further  experiments  by  others  [10,  11]  indicated 
that  these  neurons  were  not  line  detectors  but  rather 

spatial  frequency  filters.  Some  models  have  already  used 
this  fact  [12], 


The  above  discussion  motivates  the  following: 

Pef inition . The  Mandala  transform  is  an  operator 
-*  . Let  h be  a smooth  function  such 

that  h ( 0 ) = 1 then 


/v 

Mf  (x-j  »X2»w-|  jWg)  = f(x,ou)  = 


h(E,)f(^-x)exp(-i<?  ,co>)dy(e) 


where  X,atFR-  This  is  a quite  familiar  object  in  the 

engineering  circles  around  speech  processing  [13,  14]  where, 
of  course,  functions  on  the  line  instead  of  on  the  plane  are 
considered.  The  weighting  function  h plays  the  role  of  our 
window.  Indeed,  in  the  engineering  literature,  h is  called 
the  window  function  and  is  usually  chosen  to  be  a Hamming 
window  or  Hanning  window.  If  x is  fixed  then  w°  can  see 
that  the  operation  of  the  transform  is  to  first  position  f 
under  the  window  and  then  fourier  transform  the  product. 


- i if  I ' - iMaur  ' I 
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The  inversion  formula  is 


f(x)  = ( 2 tt ) j of(x,uj)exp(  i<x,w>)du(oj) 
R 


To  define  a representation  in  the  new  space  we  first  need 
some  simple  facts. 

Proposition . Let  T be  a non-singular  linear  transformation 

[J 

of  the  plane.  Then 


(fTr(x,w)  = f (Tx,(T*)-1w)  -^tT 


A 

where  f = Mf  and  fT(x)  = f(Tx). 

This  follows  by  substituting  into  the  inversion 
formula  and  using  the  adjoint  at  the  inner  product. 

ICftEftSUion  • ( f?)  ( x » to ) = f (x  + 5,03)e1<C’w> 

where  ff(x)  = f (x+£). 


Proof.  Similar. 


The  new  representation  is  now  defined  accordingly. 
Notice  that  for  a rotation  0 we  have  (0*)"1  = 0 since  0 is 
real  unitary. 


Why  will  our  new  representation  be  any  different  than 
the  old?  This  question  points  out  a number  of  interesting 
effects  that  arise  from  the  Heisenberg  uncertainty 


principle.  In  particular,  the  uncertainty  principle  has  a 
lot  to  say  about  the  "shape"  of  points  in  our  new  space 
»V* ) that  are  images  under  the  operator  M [15]. Thus, 

because  some  functions  violate  this  3hape  criterion  the  old 
space  is  embedded  within  the  new  3pace  and  can  be  decomposed 
in  a way  not  possible  before. 

The  uncertainty  principle  also  elucidates 

mathematically  the  intuitive  notion  of  texture.  In  our  new 

A 

space  an  image  is  represented  as  a function  f(x,w)  where 

2 

both  x,  ueR  . In  the  following  sense  x is  a "feature"  or 
"object"  space  while  is  a "texture"  space. 

Features  in  a given  image  can  become  texture  when  seen 
from  afar.  Similarly  when  one  examines  a texture  pattern 
that  is  suitably  magnified,  features  appear.  These  of 
course  are  very  imprecise  notions.  Furthermore,  the 
classification  of  an  image  into  texture  >r  feature  is  quite 
easy  for  extremes  but  questionable  otherwise.  Consider  a 
texture  pattern  and  begin  zooming  in  on  it.  The  perceived 

change  is  gradual.  Any  abruptness  really  stems  from  the 
language  we  use  to  describe  these  things.  The  uncertainty 
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principle  under  this  zooming  experiment  will  cause  an  image 
to  slowly  rotate  out  of  the  feature  space  and  into  the 
texture  space. 

Another  attribute  of  the  Mandala  transform  space  is 
that  it  represents  both  local  and  global  aspects  of  the 
image.  (These  terms  are  meant  in  the  engineering  sense.) 
That  a successful  visual  model  must  necessarily  represent 

both  aspects  of  an  image  and  of  the  complex  interactions 
between  them  can  be  seen  by  citing  the  Cornsweet  illusion, 
or  the  phenomenon  of  subjective  contours  [16,  17]  in  which 
strictly  local  or  texture  aspects  will  cause  striking  global 
illusions . 

It  must  be  mentioned  also  that  the  Mandala  transform 
boars  a distinct  resemblance  to  several  procedures  well 
known  to  engineers  in  the  Image  sciences  [18,  19].  It  is 
somewhat  curious  this  transform  is  not  in  common  use,  given 
its  yeoman  service  to  the  speech  community. 

6.  HARMONIC  ANALYSIS 

Now  that  the  representation  is  defined  and  the  dual  G 
is  in  hand  we  have  only  to  construct  a direct  integral 
decomposition  [20],  Instead  of  a factor  decomposition  we 
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really  need  Mautner’s  adaptation  [21].  However,  since  our 
groups  are  of  type  Lie  we  have  a method  for  finding  the 
maximal  subalgebra  of  the  commutant.  We  follow  Bargmann 
[22]  by  finding  the  center  of  the  associated  Universal 
Enveloping  algebra  and  constructing  the  derived 
representation  with  the  aid  of  the  theorems  of  Garding, 
Segal,  Nelson  and  Stinespring  [ 2 3 » 24,  25,  26].  Thus  we  are 
able  to  construct  the  "equations  or  vision"  similar  to  the 
field  equations  of  particle  physics  [27].  Solving  these 
generalized  Laplace  equations  wo  obtain  the  kernels  for  the 

desired  projection  operators. 

In  this  way  we  are  able  to  construct  a perceptually 
important  transform  such  that  a modification  in  the 
transform  space  becomes  a perceptual  correlate. 
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One  of  the  most  difficult  and  serious  problems 
associated  with  high-resolution  optical  imaging  systems  is 
the  distortion  introduced  by  atmospheric  turbulence.  The 
bulk  of  this  distortion  is  due  to  random  phase  perturbations 
caused  by  variations  of  the  refractive  index  along  the 
propagation  path.  If  images  produced*  in  the  presence  of 
atmospheric  turbulence  are  simply  recorded  and  processed  by 
conventional  methods,  the  result  is  often  a loss  of 
essential  high-spatial-frequency  detail.  This  problem  is 
serious  even  under  favorable  conditions  at  a good  site,  and 
it  is  made  worse  by  such  effects  as  wind,  solar  heating,  and 
boundary-layer  turbulence  near  a moving  object.  Much  of  the 
potential  resolving  power  of  a large  telescope  will  be  lost 
unless  the  received  optical  field  is  recorded  and  processed 
properly . 

The  removal  of  the  effects  of  atmosDheric  turbulence 
from  optical  images  is  a problem  of  long  standing.  Recently 
three  new  techniques  have  been  developed  which  could 
potentially  lead  to  significant  improvements  in  the 
resoultion  achievable  when  imaging  through  a turbulent 
atmosphere.  In  Section  2,  we  provide  some  historical 
background  and  then  discuss  these  recent  developments  in 
some  detail.  Section  3 contains  a description  and 
discussion  of  some  one -dimensional  simulations  which  we 
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performed  to  test  the  theory  and  to  provide  an  indication  of 
the  effects  of  sensor  noise  on  the  performances  of  these 
procedures.  In  Section  4,  we  summarize  our  results,  discuss 
tha  current  status  of  the  problem,  and  provide  some 
guidelines  for  further  research. 

2 . Summary  of  Previous  Work 


In  this  section,  we  provide  a brief  summary  of  previous 
work  on  the  problem  of  restoring  images  which  have  been 
blurred  by  atmospheric  turbulence.  We  make  no  attempt  to  be 
comprehensive,  but  rather  restrict  our  at  j tion  to  work 
which  is  immediately  relevant  to  "ur  problem.  For  a more 
complete  discussion  of  other  work,  and  for  a comprehensive 
bibliography,  see  [1], 

2.1.  Long-Time  Averaging 

The  problem  in  which  we  are  interested  is  characterized 
by  an  object  intensity  distribution  which  is  constant  for  a 
long  period  of  time,  during  which  the  conditions  of  the 
atmosphere  change  appreciably.  One  obvious  approach  to 
removing  the  effects  of  random  atmospheric  turbulence  is  to 
calculate  a long-time  average  oV  the  image  intensity 
distribution.  This  type  of  averaging  can  be  effected  either 
by  exposing  film  or  another  recording  medium  for  a long 
time,  or  by  averaging  many  short-time  exposures.  In  either 
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case,  averaging  has  the  effect  of  removing  the  random 
fluctuations  from  the  resulting  image.  Figure  1 shows  such 
c.n  image. 


Unfortunately,  this  long-time  averaging  also  has  the 
effect  of  removing  or  greatly  attenuating  the  high  spatial 
frequencies  from  the  image  [2].  Thus,  although  the  average 
image  is  stable,  its  spat ial -frequency  content  is  limited  by 
atmospheric  "seeing"  rather  than  by  the  imaging  system. 
Furthermore,  only  limited  improvement  can  be  achieved  by 
applying  standard  restoration  techniques  to  the  long-time 
average  image  [33 > As  an  example  of  the  seriousness  of  this 
effect,  consider  the  ?00-imh  Hale  telescope.  The  diameter 
of  the  Airy  disk  for  this  telescope  is  about  0.05  arc 
seconds,  while  a typical  seeing  limit  on  a quiet  night  is 
about  2 arc  seconds.  Thus,  if  long-time  averaging  is 
employed.  atmospheric  turbulence  degrades  the  resolving 
power  of  this  system  by  a factor  of  about  40.  If 
high-resolution  imaging  systems  are  to  reach  their  potential 
in  the  presence  of  atmospheric  turbulenie,  then,  it  is  clear 
that  alternative  methods  are  required  which  yield 
statistical  stability  while  preserving  high  spatial 
frequencies . 


Figure  1. 


Image  blurred  by  atmospheric  turbulence.  Two 
second  time  exposure  made  through  a questar 
telescope. 


2.2  . 


Interferomet rlc  Methods 


Interferometers  have  been  used  in  high-resolution 
astronomy  since  the  beginning  of  the  century.  In  1920, 
Micheison  [4]  and  Anderson  [5]  showed  that  a stellar 
interferometer  could  be  used  to  resolve  objects  which  were 
normally  beyond  the  seeing  limit  of  the  atmosphere. 
Although  the  stellar  interferometer  is  theoretically  capable 
of  measuring  both  the  amplitude  and  the  phase  of  the 
coherence  function  and  thereby  completely  recovering  the 
object  intensity  distribution,  in  practice  the  phase  cannot 
be  accurately  determined.  This  limits  the  usefulness  of  the 
stellar  interferometer  to  regular  objects  such  as  binary 
stars  whose  Fourier  transforms  have  well-defined  zeroes. 
Furthermore,  the  stellar  interferometer  is  inherently  a 
one-dimensional  scanning  device  whose  application  to 
two-dimensional  objects  would  be  cumbersome  and  expensive. 

A modified  interferometer  known  as  the  intensity 
interferometer  was  introduced  by  Hanbur y-Brown  and  Twiss  [6] 
to  improve  upon  the  performance  of  the  stp  lar 
interferometer  in  certain  situations.  This  device  is 
inherently  restricted  to  measuring  the  amplitude  of  the 
coherence  function,  and  like  the  stellar  interferometer  is 
basically  a one-dimensional  instrument.  Thus,  for  our 
purposes,  it  suffers  from  essentially  the  same  limitations 
as  does  the  stellar  interferometer.  Hence,  alternative 
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procedures  are  needed  for  those  cases  in  which  the  object 
intensity  distribution  is  asymmetric  and  in  which  complete 
two-dimensional  Information  about  the  object  is  desired. 

A number  of  modifications  to  the  basic  stellar  and 
intensity  interferometers  have  been  suggested  over  the 
years.  We  mention  two  of  these  proposed  modifications  very 
briefly.  For  a more  comprehensive  discussion,  see  [1], 

Since  the  intensity  interferometer  by  itself  cannot 
measure  phase,  it  must  be  modified  or  augmented  if  the 
complete  object  is  to  be  reconstructed.  If  a known  object 
such  as  a point  source  is  present  in  the  field  of  view  along 
with  the  unknown  object,  the  phase  of  the  unknown  object  can 
be  obtained  from  the  cross  term  resulting  from  interference 
with  the  reference  object.  This  is  entirely  analogous  to 
the  use  of  a point  source  to  retrieve  phase  in  a hologram. 
Although  this  is  an  interesting  idea,  its  application  to 
real  astronomical  problems  is  obviously  limited  because  of 
the  required  presence  of  a known  reference  object. 

An  interesting  technique  which  combines  features  of  the 
stellar  interferometer  with  those  of  the  intensity 
interferometer  has  been  suggested  by  MacPhie  [7,  8].  In 
principle,  this  device  can  measure  both  amplitude  and  phase, 
and  can  therefore  reconstruct  the  object  intensity.  Despite 
its  apparent  advantages,  however,  this  technique  has  not 
been  applied  to  real  data.  Furthermore,  like  its 
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predecessors  the  stellar  interferometer  and  the  intensity 
interferometer,  it  is  inherently  a one-dimensional  system 
whose  adaptation  to  extended  two-dirnensional  objects  would 
be  cumbersome  at  best. 

2.3.  Speckle  jntpr f^rp.tqetry. 

A partial  solution  to  the  problem  of  reconstructing 
two-dimensional  objects  was  described  by  Labeyrie  1 9 , 10], 
whose  technique  has  come  to  be  known  as  speckle 
interferometry  because  of  the  speckled  character  of  the 
images  used.  An  example  of  such  an  image  is  shown  in  Figure 
2.  In  this  technique,  a series  of  short-time  exposures  of 
the  atmospherically-degraded  object  are  taken  through  a 
narrow-band  spectral  filter.  The  exposure  time  is  short 
enough  that  the  atmosphere  is  effectively  "frozen"  during 
the  exposure,  and  successive  photographs  are  sufficiently 
separated  in  time  that  they  may  be  taken  to  be  statistically 
independent  for  the  purpose  of  averaging. 

Instead  of  simply  averaging  these  short-time  exposures 
directly,  a process  which  leads  to  a loss  of  high  spatial 
frequencies  as  we  have  seen,  speckle  interferometry  requires 
some  additional  operations  before  averaging.  First,  the 
Fourier  transform  of  each  image  is  taken,  either  optically 
or  digitally.  The  averaging  is  then  performed  on  the 
squared  magnitude  of  these  individual  Fourier  transforms. 


Figure  2 . 
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Speckle  image  of  Arcturus  made  on  the  M a y a 1 1 
four-meter  telescope  at  Kitt  Peak  by  S.  D.  Worden 
and  B.  Baxter  Delays  imposed  on  the  wavefront 
by  its  passage  through  the  atmosphere  cause 
speckles  whose  size  is  characteri sti c of  the 
diffraction  limit  of  the  telescope. 
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The  process  of  taking  the  square  of  the  magnitude 
destroys  the  information  regarding  the  phase  of  the  image. 
In  this  fact  lies  both  the  strength  and  the  weakness  of 
speckle  interferometry.  Since  it  is  the  random  phases  of 
the  high-spatial-frequency  components  which  cause  their 
mutual  cancellation  in  the  process  of  long-time  averaging, 
the  suppression  of  phase  prevents  this  cancellation  and  thus 
preserves  some  high-frequency  information  which  would  be 
lost  under  direct  long-time  averaging.  Unfortunately, 
information  about  the  phase  of  the  object  is  also 
suppressed.  This  loss  of  phase  information  limits  the  class 
of  objects  to  which  speckle  interferometry  can  be  usefully 
applied.  For  example,  this  technique  can  yield  useful 
information  about  symmetric  objects  or  binary  star  systems, 
but  not  about  an  object  with  an  arbitrary  intensity 
distribution  or  arbitrary  shape. 

Labeyrie  and  his  colleagues  have  verified 
experimentally  that  speckle  interferometry  yields 
essentially  diffraction-limited  information  about  such 
objects  as  binary  star  systems.  In  addition,  theoretical 
discussions  of  this  technique  have  been  provided  by  Miller 
and  Korff  [11]  and  by  Korff  [12].  These  authors  have  used 
models  of  propagation  in  a turbulent  atmosphere  to  explain 
why  speckle  interferometry  work3  and  to  generalize  the 
results  to  par tially-coherent  objects. 
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2.4.  The  Knox-Thompson  Technique 

The  success  of  speckle  interferometry  in  reducing  the 
effects  of  atmospheric  turbulence,  coupled  with  its  limited 
applicabilitv f has  spurred  a search  for  generalizations  and 
improvements  of  this  technique.  A significant  extension  of 
sreckle  interferometry  was  developed  by  Knox  and  Thompson  in 
1974  [1].  This  new  procedure  uses  speckle  interferometry  as 
before  to  obtain  amplitude  information.  To  this  amplitude 
information  is  added  phase  information  obtained  by  averaging 
and  processing  the  Fourier  transforms  of  the  short-time 
exposures  in  a different  way. 

In  this  section  we  present  a brief  discussion  of  the 
Knox-Thompson  technique.  This  is  not  intended  to  be  a 
complete  treatment,  but  merely  an  attempt  to  cover  the 
essential  features  of  the  procedure  in  order  to  make  the 
present  report  reasonably  self-contained.  For  a more 
thorough  treatment,  see  [1]. 
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Object  Aperture  Image 

plane  plane  plane 


Figure  3.  Imaging  configuration 
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Instantaneous  (short-exposure) 
point  spread  function  of 
atmosphere-telescope  combination. 
Fourier  transforms  of  I , I, 
and  S,  respectively. 

Focal  length  of  telescope. 


Suppose  we  take  a single  short-exposure  photograph  of 
the  object.  Then  the  intensity  distribution  in  the  image  i.i 
given  by  [13] 


I(x) 


I («) 

0 ' ' 


J -oo 


x + — a] da 


(1) 


where  z is  the  distance  of  the  object  from  the  telescope. 
Taking  the  Fourier  transform  of  both  sides  of  (1)  yields 


I(u  ) 


(2) 


The  quantity  S(u)  is  the  instantaneous  incoherent  optical 
transfer  function  of  the  atmosphere- telescope  combination. 
In  the  above,  as  well  as  succeeding  equations,  we  assume 
one-dimensional  imaging  for  notational  convenience  *nd  we 
restrict  ourselves  to  incoherent  imaging. 
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If  we  simply  average  a large  number  of  short'  exposure 
images  (say  of  duration  10  ms),  we  have  from  (1j 


< I ( x ) > = 


IQ(a)  <S^x  + j a|  >da 


(3) 


where  <S(x)>  is  the  average  of  a large  number  of 
instantaneous  point  spread  functions.  Similarly,  in  the 
frequency  domain, 


< I ( U ) > = I 


o (P)^(U) 


(4) 


where  <S(u)>  is  the  average  transfer  function.  We  assume 
that  JQ  is  constant  throughout  the  averaging  process. 

We  have  already  pointed  out  that  long-time  averaging 
(or  equivalently,  time-exposure  photography)  results  in 
appreciable  loss  of  high-spat lal-frequency  information. 
Figure  4 illustrates  this  phenomenon  in  terms  of  the 
averaged  point  spread  function  <S(x)>.  Figure  4a  shows  the 
diffract  ion -limited  point  spread  function  corresponding  to  a 
200-inch  one-dimensional  aperture.  Figure  4b  shows  a 
typical  short-exposure  point  spread  function  for  the 
atmosphere-telescope  combination.  Figure  4c  shows  a typical 
long-exposure  (or  averaged)  point  spread  function.  Since 
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this  averaged  point  spread  function  is  of  much  greater 
extent  than  the  diffraction-limited  point  spread  function, 
the  averaged  image  will  be  much  smoother  than  the 
diffraction-limited  image  and  will  thus  contain  appreciably 
less  high-frequency  detail. 

This  phenomenon  can  also  be  viewed  in  the  frequency 
domain  as  illustrated  in  Figure  5.  Figure  5a  shows  the 
optical  transfer  function  of  the  same  aperture  as  that  used 
for  Figure  4.  Figure  5b  shows  a typical  short-exposure 
transfer  function  for  the  atmosphere-telescope  combination, 
and  Figure  5c  shows  an  averaged  transfer  function. 
Comparison  of  Figure  5a  with  Figure  5c  shows  the  loss  of 
high  spatial  frequencies  incurred  oy  long-time  averaging. 
For  many  purposes,  this  is  clearly  undesirable. 

In  speckle  interferometry,  the  squared  magnitude  of  the 
image  Fourier  transform  is  calculated  before  the  averaging 
is  performed.  Performing  this  operation  on  (2)  yields 


i \mmml  i 


Page  37 


Labeyrie  has  3hown  experimentally  that<|S(u)|  > has  useful 
signal- to-noise  ratio  out  to  the  diffraction  limit  of  the 
telescope.  Korff  [12]  and  Korff  and  Miller  [11]  have 
provided  theoretical  explanations  and  have  extended  the 
results  to  partially  coherent  illumination.  A typical  plot 
of  <|S(u)J  > is  shown  in  Figure  13a.  Comparison  of  Figure 
13a  with  Figure  5c  illustrates  the  potential  gain  in 
high-frequency  information  using  speckle  interferometry. 


If  we  can  locate  a point  source  of  light  not  too  far 
removed  from  the  object  (but  not  necessarily  in  the  same 


isoplanatic  patch),  we  can  take  a second  series  of 
short-exposure  photographs  of  this  point  source,  Fourier 
transform,  square,  and  average  to  determine  < | S ( u ) | > . From 

(5)  we  will  then  have 


<|s(u)i2> 


In  practice,  we  may  wish  to  modify  this  procedure  to  avoid 
undue  amplification  of  sensor  noise,  but  the  basic  idea  is 


contained  in  ( 6 ) . 
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Labeyrie  and  hi3  colleagues  [9,  10]  have  applied 
speckle  interferometry  to  symmetrical  objects  with 
considerable  success.  Since  no  phase  information  is 
obtained  using  this  technique,  however,  it  will  not  be 
effective  if  applied  to  arbitrary  objects.  To  overcome  this 
limitation,  Knox  and  Thompson  [1]  developed  a method  foi 
calculating  the  phase  of  I (u)  by  processing  thr  same  set  of 
short-exposure  photographs  in  a somew..at  different  way.  We 
now  turn  to  a discussion  of  their  procedure. 

The  effective  pupil  function  of  the 
atmosphere-telescope  combination  is  the  product  of  the 
telescope  pupil  function  of  A(x)  and  the  instantaneous 

A 

wavefront  V(x).  In  order  to  obtain  further  results,  we  must 
make  some  assumptions  about  V(x).  It  has  been  shown  that 
amplitude  perturbations  caused  by  atmospheric  turbulence  are 
less  serious  than  the  effects  of  phase  perturbations.  Thus, 

to  a good  approximation,  we  may  neglect  the  amplitude 
variations  and  approximate  the  effects  of  the  atmosphere  by 
a phase  perturbation: 


V ( X ) 


= e 


i $ ( x ) 


(7) 


^ - 


Page  39 


We  further  assume  that  <j>(  x ) is  a Gaussian  random  process 
with  zero  mean  and  autocorrelation  function 

< 4> ( x -j  ) 4> ( x 2 ) >=  a2 k2  P ( x 2 " x-|)  f8) 

This  assumption  implies  that  we  have  removed  any  gross  phase 
shift  by  centering  the  photographs  before  processing. 

The  instantaneous  optical  transfer  function  can  now  be 
written  as  the  autocorrelation  function  of  the  effective 
atmosphere- telescope  pupil  function  C 1 3 3 • Combining  this 
result  with  the  above  assumption  yields 


Using  this  equation  and  the  assumed  properties  of  the  phase 
process  ( x ) , we  can  calculate  the  properties  of  S(u)« 


Although  we  have  alreadv  discussed  the  long-time 

w~ 


average  transfer  function  <S(u)>,  it  is  instructive  to 


calculate  this  quantity  from  (9)  using  the  properties  of 

A 

4>(x).  Averaging  both  sides  of  (9)  and  interchanging  orders 
of  expectation  and  integration  on  the  right  side  yields 
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<S(u  )> 


i <t>  ( x ) - i 4> 


<e 


(x+  { u) 


> A ( x ) 


(10) 


The  average  inside  the  integral  is  just  the  joint 
characteristic  function  of  the  Gaussian  random  process  <t>(x) 
and  is  given  by 


<e 


~ f \ 

1 4>  ( x ) - i <p [x  + ^ u) 


2.2  2,  2d 

-o  k +o  k P 


> = e 


(£ u) 


(11) 


Defining 


K(x) 


2.2  2.  2D/% 

k +o  k P ( x ) 


(12) 


and  substituting  (11)  and  (12)  into  (10)  yields 


00 

< S ( u )>  = u)  j A(x)  A (x  + £ u)  dx 


(13) 
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An  interesting  feature  of  (13)  is  the  fact  that  the 
long-time  average  transfer  function  <S(u)>  appears  as  the 
product  of  two  factors,  one  depending  only  on  the  telescope 
and  the  other  depending  only  on  the  atmosphere.  The 
integral  appearing  in  (13)  i3  just  the  transfer  function  of 
the  telescope,  while  40  represents  the  effect  of  the 
turbulent  atmosphere.  For  a high-resolution  telescope,  K 
will  be  many  times  narrower  that  the  diffraction-limited 
OTF,  and  will  thus  seriously  limit  the  usefulness  of  <S(u)>. 


The  central  feature  of  the  Knox-Thompson  procedure 
involves  the  determination  of  the  phase  of  the  object 
transform  from  the  autocorrelation  function  of  the  image 
transform.  To  see  how  this  is  done,  first  recall  (2): 


i<“>  ■ \ (r)  s(u) 


Calculating  the  autocorrelation  function  of  both  sides  of 
(2)  yields 


<I(U,)  I (u2)>  - I0 


♦ ru  i 


•<S(u-j)  S (u2)> 


--  t-  u 


where 


IQ(u)  = U0(u)l  e 


ie(u) 


Solving  ( 1 1* ) for  the  phase-difference  term  yields 


z )"ieCrO  <*(“!>  I*(“2)>  | <S(u t ) S*(u2)>| 


/fu 


T (u2)>|  <S(U| ) S ( u 2 ) > 
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(15) 


(16) 


From  (16),  we  see  that  the  difference  in  the  phases  of 


the  object  transform  evaluated  at  points  fu^/z  and  fu^ /z  can 


be  determined  from  the  autocorrelation  functions  of  the 


image  and  the  transfer  function  evaluated  at  points  u^  and 


u„  . The  phase  of  I (u)  can  be  determined  (up  to  an 
2 o 


arbitrary  constant)  by  summing  these  phase  differences. 


1 


8 

£ - 


Knox  and  Thompson  [ 1 ] have  shown  that  this  procedure  is 

effective  provided  that  the  points  u1  and  u?  are 
sufficiently  close  together  (this  statement  will  be  made 
more  precise  in  the  next  section).  Furthermore,  when  u and 
u2  are  sufficiently  close  together,  the  phase  of 
< S ( u i ) S (u2)>is  very  close  to  zero  and  can  thus  be  neglected 
without  serious  error.  The  expression  for  the  phase 
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difference  then  reduces  to 


<I(u1  ) I ( d 2 ) > 

I * ( u 2 ) > i 


(17) 


In  summary,  the  Knox-Thompson  procedure  involves  the 
following  steps: 

1.  Take  a series  of  short-exposure  photographs  of 
the  image  intensity  distribution  I(x)  through  a 
narrowband  optical  filter. 

2.  Take  a second  series  of  short-exposure 

photographs  of  a point  star. 

3.  Fourier  transform  the  photographs  in  (1)  and  (2) 
to  provide  sample  functions  of  I(u)  and  S(u), 
respectively . 

4.  Determine  the  amplitude  of  the  object  transform 
using 
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5.  Determine  the  phase  of  the  object  transform  by 
summing  the  phase  differences  obtained  from 


i6(^)'i0(-yO  <I(u])  I*(u2)> 

i (u2)>| 


6.  Take  the  inverse  Fourier  transform  of  I (u)  to 

o 

find  the  object  distribution  I (a)  . 

This  discussion  of  the  Knox-Thompson  technique  is  by  no 
means  complete.  The  reader  interested  in  more  detail  should 
consult  [1],  Some  additional  quantitative  detail  is 
provided  in  Section  3 of  the  present  report  as  we  discuss 
our  experimental  procedures  and  results. 

Knox  [1]  conducted  a series  of  one-dimensional 
simulations  of  his  procedure  in  the  absence  of  sensor  noise 
and  verified  the  essential  features  of  its  performance.  We 
have  conducted  a similar  series  of  simulations,  and  have 
included  the  effects  of  sensor  noise.  These  results  are 
described  in  Section  3 of  this  report,  our  simulations  and 


4 


if 
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those  of  others  [14]  have  shown  that  the  Knox-Thompson 
procedure  is  very  sensitive  to  sensor  noise.  These  effects 
can  be  mitigated  somewhat  by  modifying  the  procedure  to 
account  for  noise  (e.g.,  use  Wiener  filtering  to  estimate 
amplitude  and  phase),  but  the  basic  problem  remains.  A 
recently-developed  alternative  to  the  Knox-Thompson 
procedure  seems  to  promise  greater  immunity  to  sensor  noise. 
This  proceduure  is  briefly  described  in  the  next  section. 

2.5.  Sherman ' 3 Method 


The  Knox-Thompson  procedure  obtains  amplitude  and  phase 
information  separately  from  two  different  averages  taken 
over  a number  of  short-time  images.  The  amplitude  and  phase 
are  then  combined  and  the  inverse  Fourier  transform  taken  to 
obtain  an  estimate  of  the  object  intenity  distribution.  An 
alternative  procedure  which  treats  the  object  as  a whole 
without  explicitly  separating  amplitude  and  phase  was 
recently  described  by  Sherman  [14]. 


Sherman's  first  step  is  to  obtain  an  integral  equation 
which  describes  tnc  relationships  among  the  object  intensity 
distribution,  the  image  intensity  distribution,  the  combined 
effects  of  the  atmospheric  turbulence  and  the  telescope,  and 
any  imaging  noise  which  might  be  present.  The  ultimate 


objective  is  to  solve  this  integral  equation  for  the  object 
intensity  distribution  in  terms  of  measured  data  and  known 
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statistical  parameters. 

Sherman  proposes  the  following  procedure  for  obtaining 
an  approximate  solution  for  the  object.  He  first  takes  the 
Fourier  transform  of  the  images  and  then  calculates  the 
sample  covariance  matrix  of  the  image  transforms.  From  this 
quantity  he  estimates  the  Cartesian  product  of  the  object 
transform  using  standard  Wiener  filtering  together  with  some 
assumptions  about  signal  and  noise  statistics.  Once  an 
estimate  of  the  Cartesian  product  is  obtained, an  estimate  of 
the  object  itself  follows  directly  from  t.ie  fact  that  the 
only  nontrivial  eigenvector  of  the  Cartesian  product  is  the 
object  transform  itself.  If  the  estimate  of  the  Cartesian 
product  is  accurate,  the  eigenvector  associated  with  the 


» 

largest 

eigenvalue 

should  be  an 

accurate 

estimate 

of 

the 

object 

transform , 

although  a 

detailed 

analysis 

of 

the 

accuracy  has  not  been  made. 

Sherman  has  applied  this  procedure  to  simulated  data 
with  good  results.  These  simulations  also  indicate  that  his 
procedure  is  less  sensitive  to  noise  than  is  the 
Knox-Thompson  technique. 

There  are  at  least  three  features  of  the  Sherman 
procedure  as  described  in  [14]  which  should  be  investigated 
if  that  procedure  is  to  become  practical  and  useful.  First, 
and  perhaps  most  important,  is  the  large  computational 
burden  associated  with  the  procedure.  Second,  alternative 

. 


techniques  for  estimating  the  Cartesian  product  should  be 
investigated,  Third,  a more  adequate  noise  analysis  using 
appropriate  noise  models  is  needed.  With  appropriate 
modifications,  Sherman's  proc  'ure  might  well  prove  quite 
effective  in  combating  the  effects  of  atmospheric 
turbulence . 

2.6.  Real-Time  Ph/tse  Compensation 

An  entirely  different  approach  to  the  reduction  of 
atmospheric  turbulence  effects  is  to  try  to  deal  with  these 

effects  before  an  image  is  recorded  rather  than  after.  In 
many  circumstances,  the  effects  of  atmospheric  turbulence 
can  be  modeled  as  a phase  perturbation  of  the  received 
wavefront  at  the  imaging  aperture.  If  this  phase 
perturbation  were  known,  or  if  an  accurate  estimate  could  be 
obtained,  the  effects  of  turbulence  could  be  removed  by 
phase  compensation.  Several  investigators  have  considered 
the  possibility  of  implementing  real-time  phase-compensation 
systems  [15,  16].  In  these  systems,  measurements  of  optical 
phase  perturbation  are  made  over  the  objective  pupil,  either 
explicitly  or  implicitly,  with  the  result  being  usel  to 
control  a deformable  mirror  or  some  other  phase-correcting 
device  to  compensate  for  atmospheric  turbulence  in  real  time 
before  the  image  is  recorded.  Such  systems  avoid  the  heavy 
burden  of  digital  computation  which  we  have  described  above 
in  connection  with  the  Knox-Thompaon  and  Sherman  techniques, 
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but  on  the  other  hand  they  require  more  complicated  optical 
equipment  and  control  systems.  Furthermore,  the  accuracy  of 
phase  compensation  systems  is  limited  by  phase-estimation 
errors  due  to  the  finite  wavefront-sensor  signal-to-noise 
ratio  and  by  fitting  errors  due  to  the  finite  number  of 
spatial  modes  restored  by  the  wavefront  corrector  [17]. 
Thus,  we  believe  that  further  extensions  and  improvements  of 
both  phase  compensation  systems  and  post-processing 
techniques  of  the  type  we  have  discussed  above  should  be 
undertaken . 

2.7.  The  Work  of  Shapiro 

In  an  interesting  series  of  papers  [17-20],  Shapiro  has 
established  what  appears  to  be  a very  useful  framework  for 
the  analysis  and  design  of  optical  imaging  systems  which 
attempt  to  remove  the  effects  of  atmospheric  turbulence.  In 

[20] ,  he  develops  a normal-mode  decomposition  for  the 
turbulent  atmosphere  which  is  similar  tc  that  developed 
previously  for  free-space  imaging  by  Rushforth  and  Harris 

[21]  and  others.  He  then  applies  these  results  to  a study 
of  the  ultimate  performance  limits  on  imaging  through  a 
turbulent  atmosphere  [17,  19].  In  agreement  with  others,  he 
shows  that  it  is  possible  m principle  to  achieve 
diffraction-limited  imaging.  - 

Finally,  Shapiro  in  [18]  considers  the  conditions  under 
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which  diffraction-limited  imaging  may  be  achieved  even  when 
the  extent  of  the  object  is  such  that  the  imaging  is  no 
longer  isoplanatic;  i»e.,  when  the  atmospheric  point  spread 
function  depends  upon  the  position  of  the  point  source  in 
the  object  plane.  Since  all  the  work  discussed  in  Sections 
2.1  through  2.6  is  based  on  the  assumption  of  isoplanatic 
imaging,  and  since  in  many  practical  situations  the  imaging 
will  not  be  isoplanatic,  Shapiro's  results  are  of 
considerable  interest  and  importance. 


In  effect,  what  Shapiro  shows  is  that  in  many  cases  of 
interest,  the  object  can  be  broken  up  into  distinct 
isoplanatic  elements  such  that  the  contribution  of  each 
element  to  the  image  can  be  separated  from  the  others.  In 
principle,  this  enables  us  to  deal  with  each  isplanatic 
element  separately  and  then  put  the  individual  images  back 
together  to  form  the  total  image.  In  practice,  this 
procedure  may  be  formidable  in  its  complexity,  but  it  is 
nevertheless  of  interest  to  know  that  conditions  frequently 
exist  under  which  it  is  possible  in  principle. 


In  the  process  or  obtaining  his  results  on 
nonisoplanatic  imaging,  Shapiro  [18]  develops  a model  for 
the  atmosphere  which  may  be  useful  for  other  purposes.  He 
points  out  that  the  problem  of  optical  imaging  through  a 
turbulent  atmosphere  shares  many  features  with  the  problem 
of  transmitting  and  receiving  information  at  radio 
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frequencies  through  random  media  such  as  tropospheric 
scatter  channels.  In  fact,  he  shows  that  the  turbulent 
atmosphere  can  be  modeled  as  a wide-sense  stationary, 
uncorrelated  scatter  (WSSUS)  channel,  which  in  many  respects 
is  the  simplest  and  most  useful  scatter-channel  model  [22]. 
Shapiro  modifies  and  applies  the  results  of  [22]  to  show 
that  when  the  WSSUS  channel  is  underspread  (i.e.,  the 
delay-Doppier  product  is  less  than  one),  the  effects  of  the 
various  isoplanatic  elements  can  1°  separated  as  described 
above.  This  condition  frequently  occurs  in  practice. 

Thus,  we  have  an  example  in  which  a model  from 
random-channel  communication  theory  has  proven  very  useful 
in  the  analysis  of  optical  imaging  in  atmospheric 
turbulence.  In  aidition,  Shapiro  has  shown  that  this  model 
can  be  applied  to  optical  communication  systems  and  to 
speckle  interfere metry . We  believe  that  this  model,  with 
appropriate  modifications,  will  be  very  useful  in  extending 
the  work  described  in  the  report. 

? . 8 . S ummarv 

To  summarize,  we  have  described  in  this  section  several 

methods  which  have  the  potential  to  achieve  high-resolution 
imaging  in  the  presence  of  atmospheric  turbulence. 
Considerable  work  remains  before  this  potential  can  be 
realized  in  a practical  system,  however.  In  particular,  the 
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issues  of  noise  sensitivity,  computational  complexity,  and 
non i so pla na tic  imaging  must  be  investigated  more  carefully. 
Since  sensor  noise  will  be  present  in  any  practical  system, 
it  is  essential  that  we  undertstand  how  this  n^>ise  effects 
various  restoration  procedures,  and  that  we  know  how  to 
minimize  these  effects.  The  issue  of  c putational 
complexity  may  determine  whether  a technique  is  hopelessly 
impractical  or  can  be  made  operational.  Finally,  a system 
which  is  to  be  of  use  in  many  practical  situations  must  be 
able  to  deal  with  large  objects  for  which  the  imaging  will 
very  likely  be  nonisoplanatic . We  expect  the  communication 
theory  model  described  in  this  section  to  be  useful  in  these 
investigations  »s  well  as  in  suggesting  alternative  or 
modified  procedures  with  improved  performance  and 
ef f ic iency . 

3.  Simulation  of  Knox-Thompson  Procedure 

3.1.  Moiae  Simulation  with  the  Knox-Thompson  Technique 

A one-dimensional  simulation  of  the  Knox-Thompson 
procedure  for  restoring  atmospheric-turbulence-degraded 
images  was  conducted  to  obtain  a better  quantitative 
understanding  of  its  performance  in  the  presence  of  noise. 
This  work  utilized  existing  models  of  atmospheric  turbulence 
and  of  sensor  noise.  Previous  simulations  of  this  technique 
using  data  free  of  sensor  noise  have  led  to  results  in 
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general  agreement  with  predictions.  However,  it  has  been 
indicated  that  the  technique  may  be  quite  sensitive  to 
noise.  It  was  our  purpose,  therefore,  to  determine  what 
signal-to-noise  ratio  was  necessary  to  yield  an  adequate 
reconstruction  of  a degraded  image. 


3.2. 


of  the  One-Dimensional  Knox-Thompson 


This  technique  uses  the  assumption  that  the  distortions 
introduced  by  the  atmospheric  turbulence  can  be  described  by 
a Gaussian  phase  model  as  discussed  in  Section  2.4.  This 
model  has  been  used  by  several  other  authors  to  analyze  the 
effects  of  atmospheric  turbulence.  It  is  also  assumed  that 
phase  distortions  are  the  principal  mechanism  responsible 
for  the  image  degradation.  The  effect  of  the  atmosphere  was 
simulated  by  randomizing  the  phase  of  the  pupil  function  of 
the  imaging  system. 

An  array  of  Gaussian  random  noise  variables  with  unit 
variance  was  generated.  This  was  filtered  by  the  technique 
described  in  Knox  [1]  to  yield  correlated  noise.  The 
correlation  length  corresponded  to  a turbulence  cell  size  of 
12  inches  across  a 200-inch  telescope.  Figure  6 shows  128 
independent  arrays  of  random  numbers,  each  array  being  shown 
left  to  right  and  successive  arrays  top  to  bottom.  The 
arrays  were  obtained  from  a random  process  which  was  zero 
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mean,  unit  variance,  and  uncorrelated.  Tho  impulse  response 
of  the  digital  filter  is  Si  "u’n  in  Figure  7.  The  correlated 
noise  is  hown  in  Figure  8.  The  correlated  noise  arrays 
were  used  as  the  phase  (Jj^x)  and  the  aperture  function  A(x) 
shown  in  Figure  9 as  the  magnitude  to  form  the  complex  pupil 
functions.  The  pupil  functions  were  inverse  Fourier 
transformed  and  the  squared  magnitude  formed  yielding  a set 
of  128  different  point  spread  functions  S^(x),  k = 1,  2, 
. . . , 128  . 


Sk(x) 


i4>k(x)^ 


2 


The  array  of  128  point  spread  functions  is  shown  in  Figure 
10.  Each  function  is  shown  left  to  right  and  successive 
functions  are  shown  top  to  bottom. 

A second  set  of  128  point  spread  functions  was 
generated  in  a similar  manner.  These  were  convolved  with 
the  double  star  image  shown  in  Figure  11  to  yield  a set  of 
128  blurred  images  I^(x).  One  typical  blurred  image  is 
shown  in  Figure  12.  These  blurred  double  star  images  were 
used  as  the  input  for  the  Knox-Thompson  restoration. 

Let  a tilde  represent  the  Fourier  transform  of  a 


function.  Thus  S and  I represent  the  Fourier  transform  of  S 
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and  I,  respectively. 

Both  Labeyrie's  technique  and  the  Knox-Thompson 
technique  estimate  the  magnitude  of  the  Fourier  transform  of 
the  image  |IQI  in  the  same  manner.  The  first  step  is  to 
average  the  magnitude  squared  of  the  Fourier  transform  of 
the  point  spread  functions;  i.e.,  to  form  the  quantity 
<|S^(u)|  > , shown  in  Figure  13a.  The  brackets  represent 
an  average  over  the  128  spread  functions.  In  a similar 
manner,  the  average  of  the  magnitude  squared  of  the  Fourier 
transform  of  blurred  images  is  formed,  i.e.,  < | I ^ ( u ) | >. 
This  quantity  is  shown  in  Figure  13b.  Then  the  estimate  of 
the  magnitude  of  the  image  transform  |IQ|  is  given  by 


If  the  inverse  Fourier 


transform  of  j I 

o 


is  taken 


using  zero  phase,  the  autocorrelation  of  the  object 
obtained.  This  is  the  Labeyrie  result,  and  an  example 
shown  in  Figure  13c. 
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Figure  13c  Labeyrie  reconstruction  of  a double  star 
obtained  by  dividing  Figure  13b  by  Figure  13a, 
taking  the  square  root,  and  then  inverse 
transforming  using  zero  phase. 





. 


~ ZLL 


Page  61 

Knox  and  Thompson  have  gone  one  step  further  and 
provide  a method  to  obtain  phase  information.  They  have 
shown  that  differential  phase  may  be  obtained  from  the 
quantity  < I ( u ^ ) I ( U 2 ) > , where  u-j  and  u2  are  spatial 
frequencies  differing  from  each  other  by  a small  amount.  In 
our  case,  they  differed  by  one  sample  in  the  frequency 
domain . 

The  technique  is  to  average  the  differential  phase  to 
produce  an  object  phase  distribution.  This  is  combined  with 
the  amplitude  of  the  object  transform  previously  obtained 
and  inverse  transformed.  The  result  is  the  object  intensity 
distribution  with  spatial  frequencies  present  out  to  the 
diffraction  limit  of  the  telescope.  This  technique  ti.us 
offers  the  potential  to  improve  the  resolution  of  arbitrary 
intensity  distributions. 

A detailed  description  of  the  calculation  of  the  phase 
will  now  be  given.  The  first  step  is  to  shift  each  point 
spread  function  and  each  blurred  image  such  that  they  have 
no  movement  about  their  center.  Th^s  step  is  to  remove 
large  linear  phase  shifts  and  has  been  shown  [1]  to  greatly 
improve  the  contrast  of  the  resultant  reconstruction. 

The  centered  blurred  images  are  then  Fourier 
transformed.  The  transforms  are  shifted  one  step  and 
complex  conjugated.  The  product  of  the  original  and  shifted 
transforms  is  formed  and  the  average  over  the  128  images  is 
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taken.  In  mathematical  terms,  the  quantity  <I(u^)  I ( u 2 ) > 
is  obtained.  In  a similar  manner,  the  quantity 
<S(u^)  S (u2)>  is  obtained.  The  phase  of  the  quantity 

< I ( u -j  ) I *(u2)> 

< S ( u , ) S*(u2)> 


is  just  the  phase  difference  between  adjacent  points  u^  and 
u2  in  the  object  transform  array. 

The  phase  of  the  object  transform,  <jjQ  , is  obtained  by 
summing  the  phase  differences  from  the  origin  outward  (see 
Figure  14a).  This  phase  is  combined  with  the  quantity  |I  I 
and  the  inverse  Fourier  transform  taken: 


The  resultant  deblurred  image  is  shown  in  Figure  14b. 

This  result,  with  no  sensor  noise,  shows  a good 
representation  of  the  original  double  star.  The  orientation 
of  the  large  star  relative  to  the  small  star  was  reversed  as 


predicted 
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In  order  to  obtain  a more  quantitative  understanding  of 
the  effects  of  sensor  noise  on  the  Knox-Thompson  procedure, 
we  performed  a series  o.  simulations  in  which  controlled 
amounts  of  sensor  noise  were  present.  The  Knox-Thompson 
restoration  procedure  was  then  applied  to  the  noisy  data, 
and  the  resulting  restored  i radges  were  compared  with  those 
obtained  in  the  absence  of  sensor  noise. 


The  sensor-noise  model  we  used  to  generate  our  noisy 
images  is  based  on  a semiclass ical  approach  to  photon 
detection  [231.  In  this  model,  the  number  of  photoelectrons 
released  from  a small  region  of  area  A A centered  at  a point 
(x,  y)  in  the  image  plane  is  taken  to  be  a Poisson  random 
variable  with  parameter 


A = ILJlXj-JLL  t AA  (18) 

h v 

In  this  expression,  T)  is  the  quantum  efficiency  of  the 
photodetector,  h is  Planck's  constant,  v is  the  mean  optical 
frequency,  x is  the  integration  (exposure)  time,  and  I(x,  y) 
is  the  image  intensity  at  point  (x,  y). 
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It  is  well  known  that  both  the  mean  and  the  variance  of 
the  above  Poisson  random  variable  are  equal  to  A.  Thus,  the 
mean  photoelectron  current  produced  by  a detector  at  point 
(x,  y)  is  proportional  to  the  image  intensity  at  that  point, 
but  so  are  the  fluctuations  in  that  current.  If  we  define 
the  signal-to-noise  ratio  of  the  image  at  point  (x,  y)  as 
the  ratio  of  the  square  of  the  mean  current  to  the  variance 
of  the  current,  we  find  that  this  signal-to-noise  ratio  is 
just  A.  We  can  vary  A,  and  therefore  the  signal-to-noise 
ratio,  by  varying  x (the  integration  time)  or  AA  (the  area 
over  which  the  image  is  averaged). 

An  additional  simplification  results  if  we  assume  that 
the  number  of  photoelectrons  is  large.  In  this  case,  the 
photoelectric  current  will  be  approximately  Gaussian  with  a 
mean  equal  to  its  variance.  We  made  this  assumption  in  the 
simulations  which  we  performed.  The  signal-to-noise  ratio 
which  we  ascribe  to  a given  simulated  image  is  simply  the 
maximum  of  the  point-by-point  ratios  of  squared  mean  to 
variance  as  described  above.  This  definition  is  somewhat 
arbitrary,  of  course,  but  this  is  not  a serious  problem 
since  we  are  interested  primarily  in  relative  performance  as 
we  vary  the  noise. 

Results  of  the  Knox-Thompson  restoration  procedure  for 
various  signal-to-noise  ratios  as  defined  above  are  shown  in 
Figures  15a  through  1 5 j . Examples  of  Labeyrie 
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speckle-interferometry  reconstructions  are  shown  in  Figures 
16a  through  1 6 d . 


a 


Page  71 


t j 


1 • 


3.4.  Conclusions  on  the  Effects  of  Noise 


Figures  15a  through  1 5 j show  the  results  of  the 
Knox-Thompson  reconstruction  of  the  blurred  double  star 
image  with  sensor  noise.  The  signal-to-noise  ratio  (SNR) 
varied  from  5dB  up  to  70dB.  Figure  15a  has  a SNR  of  5dB. 
There  is  no  indication  of  a second  star.  Figures  15b  with  a 
10dB  SNR  could  possibly  indicate  a second  star.  Figure  15c 
through  1 5 f with  SNRs  varying  from  13dB  to  30dB  definitely 
show  a second  star,  but  the  magnitudes  of  both  the  larger 
and  smaller  stars  are  distorted.  It  is  not  until  a 40dB  SNR 
shown  in  Figure  15g  that  the  magnitudes  approach  their  true 
value.  For  higher  SNRs  (50dB  to  70dB),  there  appears  to  be 
little  improvement  in  the  reconstructed  images. 

For  SNRs  of  30dB  and  greater,  a third  "star"  appears  to 
the  right  of  the  main  star.  Its  magnitude  is  about 
one-tenth  the  magnitude  of  the  main  star.  Increasing  the 
SNR  did  not  reduce  its  presence,  but  rather  made  it  more 
distinct. 
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Figures  16a  through  1 6 d show  Labeyrie 
SNRs  varying  from  5dB  to  70dB.  Figure 
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SNR  shows  little  or  no  presence  of  a second  star.  Figure 
1 6 b with  a SNR  of  13dB  definitely  shows  a second  star,  but 
the  magnitudes  are  incorrect.  For  a SNR  of  40dB,  the 


I I 


magnitudes  are  correct. 


Increasing  the  SNR  to  70dB  shows 
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SECTION  III 

IMAGE  UNDERSTANDING 
Martin  E.  Newell 

1 . Introduction 


The  work  done  from  mid  1976  to  March  1977  in  the  area 
of  Image  Understanding  continued  to  produce  the  tools 
necessary  for  implementing  an  analysis  by  synthesis  image 
processing  facility.  The  analysis  of  images  is  carried  out 
for  the  purposes  of  automatic  recognition  of  previously 
known  objects,  or  for  synthesizing  models  of  previously 
unknown  objects.  The  basic  notion  is  that  such  analyses  can 
benefit  greatly  if  carried  out  in  conjunction  with 
three-dimensional  models  of  the  objects  in  the  scene. 

Given  modelling  and  image  synthesis  facilities  of 
sufficiently  advanced  capability,  analysis  of  such  scenes 
can  be  carried  out  using  an  analysis  by  synthesis  approach. 
Based  on  some  hypothesis  about  the  objects  in  the  scene,  and 
their  orientations  with  respect  to  the  camera,  a synthetic 
image  can  be  created  and  compared  with  the  actual  image. 
The  model  of  the  objects  in  the  scene  is  then  modified  based 
on  differences  between  these  two  images,  and  the  cycle 
repeated . 
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2.  Main  Subproblems 


Four  main  problem  areas  can  be  identified  in  attempting 


to  develop  such  a system: 


1.  Abstraction  of  perceptually  relevant  information 


from  images,  for  use  in  making  comparisons  between 


real  and  synthesized  images 


Generalized  correlation  in  both  2-D  and  3-D  for  the 


purposes  of  finding  the  best  fit  between  the  real 


image  and  a synthetic  one. 


3.  Synthesis  of  high  fidelity  images,  capable  of 


reproducing 


perceptually 


important 


characteristics  of  real  images. 


4.  Advanced  modelling  systems,  for  storing  and 


manipulating  a wide  variety  of  object 


representations. 


The  progress  that  has  been  made  in  these  four  areas  is 


described  below. 


2.1.  Abstraction  of  Perceptually  Important  Information 


In  any  Image  Understanding  task,  it  is  necessary  that 


the  features  of  an  linage  relevant  to  the  perception  of  the 


subject  of  the  image  be  readily  identified.  The  most 


obvious  and  commonly  used  features  of  images  are  high 


contrast  boundaries,  or  edges.  Several  heuristic  techniques 


for  detecting  edges  in  images  have  been  developed,  some  more 


successful  than  others  (Hueckel  (1971)[1],  Griffith 


■ 
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All  of  these  techniques  suffer  to  a greater  or  lesser 
degree  from  the  effects  of  noise,  or  of  widely  varying  local 
contrast  in  the  image.  In  an  attempt  to  provide  better 
techniques  for  abstracting  perceptually  important 
information  from  images,  a new  analysis  method  has  been 
developed.  This  new  analysis  involves  the  creation  of  a 

I 

transform  of  the  image,  called  a Mandalagram. 


A Mandalagram  is  the  two-dimensional  analog  of  a short 
time  spectogram,  used  in  the  analysis  of  acoustic  waveforms. 
A Mandalagram  is  constructed  by  fragmenting  the  given  image 
into  an  array  of  small  overlapping  areas,  typically  8x8 
picture  elements.  The  two-dimensional  Fourier  transform  of 
each  of  these  areas  is  calculated.  The  Mandalagram  is  then 
constructed  by  taking  corresponding  components  from  all  of 
the  small  area  Fourier  transforms,  and  assembling  them  into 
a mosaic  of  images,  one  image  for  each  component  of  the 
transforms,  each  component  of  the  transform  being  located  in 
its  image  in  a position  corresponding  to  the  position  of  its 
small  area  in  the  original  image. 


The  mathematics  underlying  the  Mandala  transform  is 


1 


described  in  Kajiya  (1976)[4], 
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2.2.  Generalized  Correlation 

The  ultimate  goal  of  this  part  of  the  work  is  to 
provide  correlation  techniques  to  enable  the  subject  of  a 
given  image  to  be  matched  with  one  of  a set  of  given 
prototype  objects.  The  prototype  objects  may  be 
parameterized  with  continuously  variable  parameters;  in 
which  case  the  values  of  the  parameters  to  give  a be3t  fit 
to  the  imaged  object  are  sought. 

The  problem  has  been  approached  in  several  steps  of 
increasing  difficulty.  In  the  first  step  the  given  image 
and  a prototype  image  are  assumed  to  differ  with  four 
degrees  of  freedom,  namely  x position,  y position,  rotation 
in  the  plane  of  the  image,  and  scale.  Furthermore,  the 
given  image  is  assumed  to  show  a single  object  against  a 
black  background. 

This  problem  can  be  approached  in  two  ways  - as  a 
linear  programming  problem  in  four  variables,  or  by  using 
transform  techniques  to  compute  the  correlation  function. 
Two  versions  of  the  latter  are  used  here. 

By  way  of  review,  two  images  which  differ  in  only  x 
position  and  y position  can  be  cross  correlated  using  the 
two-dimensional  Fourier  transform.  Let  the  given  and 
prototype  images  be  f(x,y)  and  g(x,y)  respectively.  Then 
the  correlation  of  the  two  can  be  computed  as  the 
two-dimensional  convolution: 


j 


C(x,y)  = f (x,y)  * g(-x,-y)  = <f,g> 
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The  juxtaposition  of  the  two  images  giving  a best  fit  is 
given  by  the  coordinates,  (x,y)  for  which  c(x,y)  has  a 
maximum.  Morever,  the  magnitude  of  the  maximum  gives  an 

indication  of  the  closeness  of  the  fit,  unity  indicating  a 
perfect  fit. 

Whereas  higher  dimensional  versions  of  the  above 
technique  can,  in  principle,  be  applied  in  the  four  degrees 
of  freedom  case,  the  required  computational  time  and  space 
become  wholly  excessive.  Consequently,  two  other  approaches 
which  avoid  such  problems  have  been  developed.  The  first  of 
these,  called  the  Transform  Method,  transforms  the  images 
into  a space  where  the  effects  of  translation  can  be 
separated  from  the  effects  of  scale  and  rotation.  This 
allows  the  problem  to  be  considered  as  two  two-dimensional 
correlations  in  se  ies.  The  secoud  alternative  approach 
called  the  Centroid  Method,  is  similar  in  intent,  is 
computationally  more  efficient,  but  more  prone  to 
inaccuracies  in  the  presence  of  noise. 
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The  Transform  Method 


Consider  the  Fourier  transforms  of  two  images  f(x,y) 


and  g(x,y) 


F(w-(,w2)  = F(f(x,y)) 


G(W|  ,w,?)  = F(g(x,y) ) 


Assume  that  f and  g differ  by  a translation  of  (Ax, Ay),  a 

rotation  in  the  x,y  plane  by  a,  and  a scale  s,  i.e. 


g(.<,y)  = f(trans(Ax,Ay)rotate(a)scale(s)(x,y)) 


then  it  may  be  shown  that: 


G(w, ,w 0)  = — F(rotate)(a)scale(l/s)(x,y))e^wlAx+w2Ay^ 

1 2 hi 


From  this  expression  we  see  that  the  transform  is  rotated 
exactly  the  same  as  the  image,  it  is  scaled  by  the  inverse 
of  the  image  scale,  and  its  phase  is  altered  by  the  addition 
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of  planar  phase.  That  is,  all  three  geometric 
transformations  affect  phase,  but  only  rotation  and  scale 
affect  magnitude,  in  the  manner  shown.  Consequently,  the 
differences  in  the  magnitudes  of  the  two  Fourier  transforms 
are  entirely  due  to  rotation  and  scale.  Moreover,  the 
magnitudes  will  themselves  differ  by  the  same  rotation  and 
inverse  scale  as  do  the  images. 

By  working  with  the  magnitudes  of  the  Fourier 
transforms,  the  problem  is  reduced  to  one  having  only  two 
degrees  of  freedom,  rotation  and  scale.  This  can  be  solved 
as  a two-dimensional  correlation.  However,  it  is  first 
necessary  to  change  the  coordinate  system  to  radial 
coordinates,  (r,0)  with  exponential  radial  distances,  i.e. 


|F(w-|  ,w2)  I -*■  <J>(r»0) 
|G(w-j  ,w2)  | -*•  y(r,6) 


Consider  the  effect  of  scaling  by  1/s  in  w*|  and  w 


|F(w1/s,w2/s)  | ■+  4>(rs,<j>s) 


where 


= In  J(w1/s)2  + (w2/s) 
= In  7w^  + w2)  (1/s)2 

Vi 


rs  = 1 


~T~  ; 

n w,  + w, 

I t 


-Ins 


= r - 1 n s 


-i  Vi 

w2/s 


-1  W1 
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i.e.  A scale  by  (1/s)  in  (Wj.Wg)  space  is  mapped  into  a 
translation  by  -In  s in  the  r coordinate  of  (r,9)  space.  It 
can  be  shown  that  a rotation  by  a in  (w-|,W2)  space  maps  into 
a translation  by  a in  the  9 coordinate  of  (r,9)  space.  This 
enables  us  to  again  use  two-dimensional  Fourier  transform 
techniques  to  cross  correlate  the  two  functions.  However, 
the  functions  4>  ( r , 9 ) and  y(r,9)  must  first  be  compensated 
for  the  distortion  of  area  that  occurs  in  the  coordinate 
transf ormation . This  step  is  necessary  because  the 
correlation  is  essentially  an  integral  with  respect  to  r and 
9 and  is  therefore  affected  by  changes  of  metric.  This 


compensation  can  be  achieved  by  multiplying  both  4>  ( r , 9 ) and 

2 


Y ( r 9 ) by  e which  is  the  determinant  of  the  Jacobian  of 
the  transformation,  i.e. 


<f> ' ( r , 6 ) = er  (Mr, 9) 


Y* (r  ,9)  = e y(r »6) 


rio  a i c o u 1 b Cl  uu  1 o uvi  i cxai/i-Ou  ^ W i lii  U V d i v v w Ui  i 


and  9 which  give  a maximum  in  the  function: 


<d'  ,j’> 


<d  r,  d ' >< j ' , j ' > 


uaiui 


^ — 


Let  these  values  be  R and  0.  Then  the  scale  factor  giving  a 
best  fit  in  the  original  image  is 

Scale  = e " 

and  the  rotation  is  0. 

Having  determined  the  rotation  and  scale,  it  remains  to 
find  the  displacements  in  x and  y.  This  can  be  done  by 
applying  the  scale  and  rotation  found  above  to  the  prototype 
image,  then  carrying  out  a two-dimensional  cross  correlation 
with  the  given  image  to  find  the  displacements.  A 
one-dimensional  version  of  this  method  is  shown  in  Figure  1. 

The  transform  method  for  determining  the  four  degrees 
of  freedom  may  be  used  in  the  presence  of  limited  amplitude 
high  frequency  noise.  Using  this  method  the  goal  of 
correlating  images  with  templates  in  2-D  can  be  achieved. 

The  Centroid  Method 

For  the  given  subDroblem  of  determining  the  position, 
size,  and  orientation  of  a given  image  against  a black 
background,  another  transform  technique  exists.  This 
involves  finding  the  centroid,  C ^ , of  the  given  image  and  Cy 
of  the  template.  The  position  of  the  centroid  is 
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independent  of  rotation  and  scale  about  the  centroid. 

Therefore  the  difference  in  the  positions  of  the  centroids 
Cj  and  C-j.  will  be  precisely  the  translation  (Ax  Ay)  between 
the  given  image  and  the  template. 

To  find  the  rotation  and  scale,  w translate  both  tho 
image  and  template  so  as  to  bring  their  centroids  to  the 
origin  of  coordinates.  From  here  the  problem  is  the  same  as 
in  the  Transform  Method,  except  that  now  we  are  in  image 
space  instead  of  in  magnitude  Fourier  transform  space. 

Since  the  computation  of  centroid  is  extremely  simple, 
the  Centroid  Method  is  faster  than  the  Transform  Method. 

However,  it  is  probably  more  prone  to  error  in  the  presence 
of  noise  in  the  image  (Figure  2). 

2.3.  High  Fidelity  Synthetic  Images 

The  basic  notion  underlying  the  present  approach  to 
Image  Understanding  is  that  it  should  be  carried  out  in  the 
context  of  a model  of  the  objects  of  interest  in  the  scene. 

That  model  will  be  created  and  refined  based  on  a comparison 
between  the  given  image  and  a synthetic  image  of  the  model. 

In  oraer  tnat  this  comparison  snouia  yieia  tne  most 
information  it  is  necessary  that  the  synthesized  images 
should  be  as  realistic  as  possible,  especially  in  those 
characteristics  that  are  perceptually  important.  Such 
characteristics  include  shape,  orientation,  lighting, 

y 
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highlights,  reflections,  surface  properties,  texture  color, 
and  shadows. 

While  current  state-of-the-art  synthetic  images  bear  a 
remarkable  resemblance  to  images  of  real  scenes,  the 
simulation  of  reflections,  surface  properties,  texture  and 
shadows  leaves  much  room  for  improvement.  To  this  end  the 
development  of  improved  capabilities  in  these  areas  has  been 
pursued  . 

Texture 

Most  image  synthesis  techniques  model  surfaces  as  being 
smooth  and  uniformly  colored.  Even  with  the  use  of 
highlighting  techniques  the  surfaces  appear  to  be  made  of 
some  plastic  type  material.  Catmull's  work  (1974)[5] 
introduced  a technique  for  mapping  patterns  onto  a surface 
for  the  purposes  of  display.  The  technique  simulates  the 
painting  of  a pattern  onto  a surface.  No  restriction  was 
made  on  the  form  of  the  pattern.  Consequently,  textured 
surfaces  could  be  simulated  by  the  use  of  appropriate 
patterns . 

As  in  all  discrete  implementations  of  continuous 
signals,  an  aliasing  problem  arises  in  the  process  of 
mapping  a pattern  onto  a surface.  This  is  a compound 
problem  in  this  case  since  the  texture  pattern  is  sampled  to 
determine  the  shade  of  a particular  piece  of  surface,  and 
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the  surface  .■*  sampled  for  the  purposes  of  display.  As  with 
all  aliasing  problems,  the  solution  is  to  band  limit  the 
signal  being  sampled  by  low  pass  filtering.  However,  the 
problem  of  band  limiting  the  mapped  pattern  is  severe  since 
it  is  sampled  at  non-uniform  intervals.  A paper  by  Blinn 
and  Newell  (1976)[6]  presents  a new  technique  for  carrying 
out  the  pattern  mapping  process  with  a minimum  of  aliasing 
problems.  This  new  technique  is  computationally  more 
efficient  than  previously  used  methods,  though  techniques 
for  carrying  out  this  type  of  process  in  real  time  have  not 
yet  been  developed.  This  paper  is  included. 


Reflections 

The  correct  simulation  of  mirror  reflections  is 
necessary  when  dealing  with  highly  polished  or  glazed 
objects.  Such  reflections  give  rise  to  highlights,  which 
are  perceptually  important  shape  clues.  Previous  attempts 
to  simulate  reflections  have  been  limited  to  the  reflection 
of  point  light  sources  in  surfaces  which  are  polished,  but 
not  mirror  reflecting. 

A new  technique  for  simulating  true  mirror  reflections 
in  curved  surfaces  has  been  developed.  This  technique  is 
also  described  in  Blinn  (1976).  Basically  the  algorithm 
functions  by  recursively  fragmenting  the  surface  in  the  same 
manner  as  is  used  in  the  patterm  mapping  algorithm.  When  a 
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fragment  is  to  be  displayed,  its  surface  normal  is  used  to 
address  a cylindrical  projection  of  the  environment 
surrounding  the  object.  The  intensity  value  so  obtained  is 
used  to  scale  the  white  specular  reflection  component  of  the 
surface  shading.  The  diffuse  reflection  component  is 
computed  using  one  of  several  previously  reported  shading 
rules,  such  as  Lambert  or  Bui-Tuong  ( 1975)C7]  . 

2.4.  Advanced  Modelling  systems 

The  development  of  an  analysis  by  synthesis  system  for 
Image  Understanding  relies  heavily  on  the  existence  of  a 
three-dimensional  model  of  the  scene  being  analyzed.  This 
model  will  be  updated  and  refined  as  the  analysis  proceeds. 

In  order  to  reduce  the  number  of  degrees  of  freedom  in 
such  a model,  the  use  of  a .single-primitive  modelling  system 
is  untenable.  For  example,  a commonly  used  orimitive  for 
modelling  three-dimensional  scenes  is  the  planar  polygon. 
However,  the  modelling  of  even  comparatively  simple  scenes 
involves  the  use  of  many  hundreds,  or  even  thousands,  of 
polygons.  Clearly,  higher  level,  parameterized,  models  are 
needed . 


The  use  of  higher  level  parameterized  models  implies 
specialization  of  the  mo  jelling  system  for  each  type  of 
model.  This  is  because  the  types  of  parameters,  and  the 
type  of  control  that  they  give  over  the  model,  will  vary 
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widely  from  one  model  to  another.  However,  such 
specialization  is  undesirable  in  that  many  different 
modelling  systems  will  be  required,  and  the  modelling  of  a 
scene  containing  more  than  one  type  of  model  becomes 
impossible . 

Consequently,  system  structures  have  been  developed 
which  permit  several  high  level  parameterized  models  to 
coexist  within  the  same  system,  and  to  communicate  with  one 
another . 

The  system  structure  that  has  been  developed  and 
implemented  is  based  on  ideas  from  several  sources, 
including  Winograd  ( 1973H8],  Hewitt  ( 1973)t9],  SIMULA 
(1973)[10],  Newell  (1975)[11].  The  essential  idea  is  that 
objects  should  be  represented  as  combinations  of  data  and 
procedure,  where  both  the  data  and  procedure  parts  can  vary 
from  object  to  object.  This  is  in  contrast  to  the  more 
conventional  use  of  only  data  to  represent  objects,  all 
objects  being  interpreted  by  a common  centralized  set  of 
procedures . 

The  added  flexibility  afforded  by  this  approach  is  just 
what  is  needed  in  a geometric  modelling  system  for  use  in 
Image  Understanding.  However,  the  implementation  of  such  a 
system  raises  many  questions,  such  as:  how  should  these 
models  be  created  and  manipulated,  at  which  level  should 
such  models  communicate,  what  facilities  in  the  system  are 
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object  data  can  be  thought  of  as  a set  of  parameters  to  the 
class  procedure,  and,  together  with  the  class  procedure  it 
defines  the  object.  The  combination  of  class  procedure  plus 
object  parameters  can  be  thought  of  as  a single-primitive 
subsystem  residing  within  the  total  system. 

Since  the  class  procedure  is  shared  among  the  objects 
of  the  class,  the  only  unique  manifestations  of  the  objects 
are  their  data  blocks.  Consequently,  the  central 
communication  and  control  procedures  access  objects  via 
their  data  blocks,  which  ii.  turn,  refer  to  the  relevant 
class  procedure.  However,  for  the  purposes  of  creating  a 
new  instance  of  a class,  the  central  procedures  must 
communicate  directly  with  the  class  procedure. 


The  nature  of  the  messages  that  are  exchanged  between 
class  procedures  and  the  central  control  procedures  needs 
careful  attention.  In  a conventional,  single-primitive 
system,  the  communication  between  the  central  procedures  and 
the  data  base  is  in  term?  of  the  primitives  of  the  system. 
However,  in  the  structure  used  here  the  communication  must 
be  in  terms  of  much  higher  level  abstractions.  For  examp.  e, 
if  it  is  necessary  to  synthesize  an  image  of  an  object  for 
the  purposes  of  comparison  with  a given  photograph,  the 
message  "give  me  your  image"  would  be  sent  to  the  object 
(class  procedure  and  data).  It  would  then  be  the  task  of 
the  class  procedure  to  create  the  image,  possibly  using  some 
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commonly  available  image  synthesis  utility  procedure,  and 
transmit  it  back  to  the  central  procedures.  The  point  is 
that  the  technique  used  to  create  the  image  is  completely 
within  the  control  of  the  class  procedure,  and  can  exploit 
any  special  information  known  about  the  class  of  objects 
represented.  Examples  of  such  information  are  symmetry, 
convexity,  and  separability. 

Another  form  of  communication  takes  place  in  the  system 
when  one  class  procedure  needs  to  access  another.  Many 
objects  in  the  real  world  are  structured  assemblies  if  other 
objects.  In  such  cases  the  class  procedure  would  make 
reference  to  the  sub-objects  by  the  same  mechanism  as  that 
• used  by  the  central  procedures.  For  example,  a request  for 

the  weight  of  such  an  object  would  result  in  a depth  first 
traverse  of  the  tree  of  sub-objects. 

The  creation  and  modification  of  objects  is  carried  out 
by  the  relevant  class  procedure.  The  creation  and 
modification  of  class  procedures  is  not  so  straightforward. 
Not  many  programming  systems  provide  for  dynamic 
modification  of  program.  Those  that  do  still  lea”j  much  to 
be  desired  at  the  human  engineering  level.  Consequently, 
the  preliminary  design  of  a new  programming  language  has 
been  completed.  Programs  written  in  this  language  would  be 
run  i t e r pr e t i ve 1 y , to  achieve  the  dynamic  capabilities. 
Moreover,  the  primitives  of  this  language  include  the 
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abstractions  dealt  with  in  the  system.  These  include 
objects,  class  procedures,  instance  blocks,  tree  structures 
of  objects,  etc.  in  addition  to  the  conventional  arithmetic 
and  control  capabilities.  This  work  is  the  subject  of  an 
ongoing  study. 

3.  Further  Work 

The  purpose  of  this  year's  work  was  to  develop  the 
tools  necessary  in  an  automated  or  semi- automated  image 
understanding  system.  While  significant  progress  was  made 
in  the  four  areas  described,  further  work  of  this  type  is 
still  needed.  In  the  area  of  abstracting  perceptually 
relevant  information,  the  properties  and  application  of  the 
Mandala  transform  need  to  be  investigated. 

In  the  area  of  generalized  correlation  the  techniques 
developed  need  extending  to  handle  multiple  objects  and 
objects  on  textured  backgrounds.  Extensions  to  the  3-D 
problem  are  also  needed.  In  this,  two  extra  rotational 
degrees  of  freedom  are  introduced.  This  invokes  the  need  to 
use  a 3-D  template  which  is  a model  of  the  object  in 
question.  For  any  trial  orientation  of  the  model  a 
synthetic  image  would  be  generated  and  correlated  with  the 
given  image  using  the  methods  described  here. 

This  brings  up  the  need  for  high  fidelity  synthetic 
images.  Again,  while  significant  advances  have  been  made, 
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especially  in  the  area  of  realistic  surface  properites, 
further  work  is  needed.  However,  the  quality  of  images  now 
being  produced  is  good  enough  for  preliminary  use  in  an 
integrated  image  understanding  system.  The  main  outstanding 
deficiency  is  still  the  simulation  of  shadows,  which  give 
important  perceptual  clues  and  should  therefore  be  included. 

The  area  of  modelling  systems  needs  little  further 
development.  The  system  that  was  implemented  was  a trial 
testbed  and  as  such  would  not  be  suitable  for  a working 
integrated  system. 

The  next  step  in  this  work  would  be  to  integrate  the 
various  parts  in  an  image  understanding  system.  The 

modelling  system  would  store  and  manipulate  3-D  models  of 
. objects  of  interest.  These  would  be  used  to  synthesize  2-D 

images  for  comparison  with  the  given  image.  Image 

enhancement  techniques  would  be  applied  to  both  the  given 
and  synthesized  images  to  concentrate  attention  on 

perceptually  important  features.  Generalized  correlation 
techniques  could  then  be  used  to  find  the  best  fit  between 
real  and  synthetic  images,  and  to  measure  the  quality  of  the 
fit.  Feedback  from  the  correlation,  either  manual  or 
automatic,  would  be  used  to  modify  the  model  and  complete 
the  iterative  cycle. 
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Introduction 

In  1974  Edwin  Catmull  [2]  developed  an  algorithm 
for  rendering  continuous  tone  images  of  objects  modeled 
with  bivariate  parametric  surface  patches.  Unlike  most 
earlier  algorithms  [6,  8 9,  10],  which  require  th  it  ob- 
jects be  approximated  by  collections  of  planar  poly- 
gons, Catmull’s  algorithm*.  works  directly  from  the 
mathematical  definition  of  the  surface  patches.  The 
algorithm  functions  by  recursively  subdividing  each 
patch  into  smaller  patches  until  the  image  of  each 
fragment  covers  only  one  picture  element.  At  this 
stage,  visibility  and  intensity  calculations  are  performed 
for  that  picture  element.  Since  the  subdivision  process 
will  generate  picture  elements  in  a somewhat  scattered 
fashion,  the  image  must  be  built  in  a memory  called  a 
depth  buffer  or  Z-buffer.  This  is  a large,  random  access 
memory  which,  for  each  picture  element,  stores  the  in- 
tensity of  the  image  and  the  depth  of  the  surface  visible 
at  that  element.  As  each  patch  fragment  is  generated, 
its  depth  is  compared  with  that  of  the  fragment  cur 
rently  occupying  the  relevant  picture  element.  If  greater, 
the  new  fragment  is  ignored,  otherwise  the  picture 
element  is  updated. 

This  paper  describes  extensions  of  CatmulTs  al- 
gorithm in  the  areas  of  texture  and  reflection.  The 
developments  make  use  of  digital  signal  processing 
theory  and  curved  surface  mathematics  to  improve 
image  quality. 


Texture  Mapping 

Catmull  recognized  the  capability  of  his  algorithm 
for  simulating  variously  textured  surfaces.  Since  the 
bivariate  patch  used  is  a mapping  of  the  unit  square  in 
the  parameter  space,  the  coordinates  of  the  square 
can  be  used  as  a curvilinear  coordinate  system  for  the 
patch.  It  is  a simple  matter  for  the  subdivision  process 
to  keep  track  of  the  parameter  limits  of  each  patch 
fragment,  thereby  yielding  the  parameter  values  at 
each  picture  element.  These  parameter  values  may  then 
be  used  as  a key  for  mapping  patterns  onto  the  sur- 
face. As  each  picture  element  is  generated,  the  para- 
metric values  of  the  patch  within  that  picture  element 
are  used  as  input  to  a pattern  definition  function.  The 
value  of  this  function  then  scales  the  intensity  of  that 
picture  element.  By  suitably  defining  the  pattern  func- 
tion, various  surface  textures  can  be  simulated. 

As  Catmull  pointed  out,  simply  sampling  the  tex- 
ture pattern  at  the  center  of  e^ch  picture  element  is 
not  sufficient  to  generate  the  desired  picture,  since  two 
adjacent  picture  elements  in  the  image  can  correspond 
to  two  widely  separated  points  in  the  patch  parameter 
space,  and  hence  to  widely  separated  locations  in  the 
texture  pattern.  Intermediate  regions,  which  should 
somehow  influence  the  intensity  pattern,  would  be 
skipped  over  entirely.  This  is  a special  case  of  a phe- 
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nomenon  known  as  “aliasing”  in  the  theory  of  digital 
signal  processing.  This  theory  [7]  treats  the  image  as  a 
continuous  signal  which  is  sampled  at  intervals  cor- 
responding to  the  distance  between  picture  elements. 
The  well-known  “sampling  theorem”  states  that  the 
sampled  picture  cannot  represent  spatial  frequencies 
greater  than  1 cycle/2  picture  elements.  “Aliasing” 
refers  to  the  result  of  sampling  a signal  containing  fre- 
quencies higher  than  this  limit.  The  high  spatial  fre- 
quencies (as  occur  in  fine  detail  or  sharp  edges)  re- 
appear under  the  alias  of  low  spatial  frequencies,  i his 
problem  is  most  familiar  as  staircase  edges  or  “jag- 
gies.”  In  the  process  of  texture  mapping,  the  aliasing 
can  be  extreme,  owing  to  the  potentially  low  sampling 
rate  across  the  texture  pattern. 

To  alleviate  this  problem  we  must  filter  out  the  high 
spatial  frequency  components  of  the  image  (in  this 
case  the  texture  pattern)  before  sampling.  This  filter- 
ing has  the  effect  of  applying  a controlled  blur  to  the 
pattern.  This  can  be  implemented  by  taking  a weighted 
average  of  values  in  the  pattern  immediately  surround- 
ing the  sampled  point.  Digital  image  processing  theory 
provides  a quantitative  measure  of  the  effectiveness  of 
such  weighting  functions  in  terms  of  how  well  they 
attenuate  high  frequencies  and  leave  low  frequencies 
intact, 

Catmull  achieved  the  effect  of  filtering  by  main- 
taining an  additional  floating-point  word  for  each 
picture  element.  This  word  contained  the  fraction  of 
the  picture  element  covered  by  patch  fragments.  For 
each  new  fragment  added  to  the  picture  element,  the 
texture  pattern  was  sampled  and  the  intensity  was 
averaged  proportionally  to  the  amount  of  the  picture 
element  covered  by  the  patch  fragment.  Examination 
of  the  spatial  frequency  filter  effectively  implemented 
by  this  technique  shows  that  it  is  much  better  than 
point  sampling  but  is  not  optimal. 

The  method  discussed  here  does  not  require  the 
extra  storage  and  uses  a better  anti-aliasing  filter. 
This  filter  is  implemented  by  a weighting  function  origi- 
nally used  by  Crow  [3]  to  minimize  aliasing  at  polygon 
edges  (“jaggies”).  It  takes  the  form  of  a square  pyra- 
mid with  a base  width  of  2X2  picture  elements,  in 


Fig.  1 - Region  of  lexture  pattern  corresponding  to  picture  element : 
left-hand  side  shows  texture;  right-hand  side  shows  image. 


the  texture  mapping  case,  the  2X2  region  surrounding 
the  given  picture  element  is  inverse  mapped  to  the  cor- 
responding quadrilateral  in  the  uy  parameter  space 
(which  is  the  same  as  the  texture  pattern  space),  see 
Figure  1.  The  values  in  the  texture  pattern  within  the 
quadrilateral  are  weighted  by  a pyramid  distorted  to 
fit  the  quadrilateral  and  summed. 

The  derivation  of  the  quadrilateral  on  the  texture 
pattern  makes  use  of  an  approximation  that  the  para- 
metric lines  within  one  picture  element  are  linear  and 
equally  spaced.  The  X,Y  position  within  a picture 
element  can  then  be  related  to  the  uy  parameters  on 
the  patch  by  a simple  affine  transformation.  This  trans- 
formation is  constructed  from  the  uy  and  X,Y  values 
which  are  known  exactly  at  the  four  corners  of  the 
patch  fragment. 

Given  this  algorithm,  we  now  investigate  the  effects 
of  various  texture  definition  methods.  We  will  use,  as 
our  sample  object,  a plain  teapot  constructed  of  26 
bicubic  patches.  First,  the  pattern  may  be  some  simple 
function  of  the  uy  parametric  values.  A useful  ex- 
ample of  this  is  a simple  gridwork  of  lines.  The  result  is 
as  though  parametric  lines  of  the  component  patches 
are  painted  on  the  surface,  Figure  2.  Note  that  the 
edges  of  the  pattern  lines  show  very  little  evidence  of 
aliasing  in  the  form  of  staircases.  Second,  the  pattern 
may  come  from  a digitized  hand  drawn  picture,  Figure 
3.  Third,  the  pattern  may  come  from  a scanned-in 
photograph  of  a real  scene,  as  in  Figure  4.  Incidentally, 
this  picture  makes  the  individual  patches  very  clear. 
This  type  of  pattern  definition  enables  the  computer 
production  of  “anamorphic”  pictures.  These  are  pic- 
tures which  are  distorted  in  such  a way  that  when  viewed 
in  a curved  mirror  the  original  picture  is  regenerated, 
Figure  3.  The  patch  itself  is  defined  so  that  the  para- 
metric lines  are  stretched  in  approximately  the  correct 
fashion  and  a real  photograph  is  mapped  onto  the 
patch.  Figure  5 should  be  viewed  in  a cylindrical  mirror 
(e.g.  a metal  pen  cap)  with  the  axis  perpendicular  to 
the  page.  The  fourth  source  of  texture  patterns  shown 
here  is  Fourier  synthesis.  A two-dimensional  frequency 
spectrum  is  specified  and  the  inverse  Fourier  transform 
generates  the  texture  pattern.  This  is  a simple  way  of 
generating  wavy  or  bumpy  patterns.  Certain  restrictions 
on  the  form  of  the  input  spectrum  must  be  followed  to 
ensure  that  the  pattern  has  an  even  distribution  of 
intensities  and  is  continuous  across  the  boundaries. 
An  example  of  this  type  of  texture  is  shown  in  Figure 
6.  T1  e texture  patterns  used  here  were  generated  before 
picture  synthesis  began  and  stored  as  256X256 
element  pictures  in  an  array  in  random  access 
memory. 


Reflection  in  Curved  Surfaces 

The  second  topic  discussed  in  this  paper  concerns 
lighting  models.  Typically,  visible  surface  algorithms 
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Fig.  2 Simple  gridwork  texture  pattern;  left-hand  side  shows  texture  pattern;  right-hand  side  shows  textured  object 


Fig  3 Hand  sketched  texture  pattern  left-hand  side  shows  texture  pattern,  right-hand  side  shows  textured  object 


Fig.  4.  Photographic  texture  pattern:  left-hand  side  shows  texture  pattern,  right-hand  side  shows  textured  object 


have  been  used  to  give  the  impression  of  shiny  surfaces, 
but  there  is  little  physical  justification  for  such  func- 
tions, and  the  range  of  effects  is  limited  The  modeling 
of  more  realistic  lighting  was  first  investigated  by  Bui- 
Tuong  Phong  [1].  His  model  of  reflection  incorporated 
a term  which  produced  a highlight  over  portions  of  the 
surface  where  the  normal  falls  midway  between  the  light 


determine  intensities  within  an  image  by  using  Lam 
bert’s  (cosine)  law-  i - s(L-N)}  where  i = intensity, 
.9  = surface  shade,  L = light  direction  vector,  N = 
surface  normal  vector,  and  denotes  vector  inner 
product. 

Variants  on  this  function,  such  as 
/ - s(L-N)**nf  for  n > 1 
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source  direction  and  the  viewing  direction.  This  is  moti- 
vated by  the  fact  that  real  surfaces  tend  to  reflect  more 
light  in  a direction  which  forms  equal  angles  of  inci- 
dence and  reflectance  with  the  surface  normal.  This  can 
be  easily  implemented  by  simulating  a virtual  light 
source  in  a direction  halfway  between  the  light  source 
and  viewing  directions,  and  raising  the  result  to  some 
high  power  to  make  the  highlights  more  distinct: 

/ - s(L-N)  + g(Lf-N)**n 

where  L!  - virtual  light  source  direction,  g = glossiness 
of  surface  (0  to  1).  Figure  7 shows  an  image  generated 
using  the  above  function  with  n = 60.  These  techniques 
work  well  for  satin  type  surfaces  but  images  of  highly 
polished  surfaces  still  lack  realism.  This  is  largely  due 
to  the  absence  of  true  reflections  of  surrounding  objects 
and  distributed  light  sources. 

The  simulation  of  reflections  in  curved  surfaces 
requires  an  accurate  model  of  the  properties  of  the 
surface  and  access  to  accurate  normal  vectors  at  all 
points  on  the  surface.  The  approximation  of  curved 
surfaces  by  collections  of  planar  polygons  is  inade- 
quate for  this  purpose,  so  extensions  of  the  techniques 
of  Gouraud  [5]  and  Bui-Tuong  Phong  [lj  hold  little 
promise. 

The  subdivision  algorithm,  however,  provides  ac- 
curate information  about  surface  position  and  can 
be  made  to  give  accurate  surface  normals  at  every 
picture  element.  This  is  the  first  algorithm  that  provides 
me  appropriate  information  for  the  simulation  of  mirror 
reflections  from  curved  surfaces.  For  each  picture  ele- 
ment, the  vector  from  the  object  to  the  observer  and 
the  normal  vector  to  the  surface  are  combined  to  de- 
termine what  part  of  the  environment  is  reflected  in 
that  surface  neighborhood.  It  can  be  shown  that,  for 
surface  normal  vector  ( Xn , Yn , Zri)  and  viewing  posi- 
tion (1,0,  0),  the  direction  reflected,  ( Xr , Kr,  Zr),  is 

Xr  = 2*Xn*Zn,  Yr  « 2 *Yn*Zn,  Zr  = 2 *ZmZn  - 1, 

Having  established  the  direction  of  the  ray  which  is 
reflected  to  the  eye,  it  remains  to  find  what  part  of  the 
environment  generated  that  ray.  For  this,  a model  of 
the  environment  is  needed  which  represents  surround- 
ing objects  and  light  sources.  Clearly,  the  view  of  the 
environment  as  seen  from  different  points  on  the  sur 
face  will  vary.  However,  if  it  is  assumed  that  the  en- 
vironment is  composed  of  objects  and  HJ  i sources 
which  are  greatly  distant  from  the  object  being  drawn, 
and  that  occlusions  of  the  environment  by  parts  of  the 
object  itself  are  ignored,  then  the  environment  can  be 
modeled  as  a nvo-dimensional  projection  surrounding 
the  drawn  object.  Stated  another  way,  the  object  is 
positioned  at  the  center  of  a large  sphere  on  the  inside 
of  which  a picture  of  the  environment  has  been  painted. 
These  simplifications  allow  the  environment  to  be 
modeled  as  a two-dimensional  intensity  map  indexed 
by  the  polar  coordinate  angles  of  the  ray  reflected 


Fig.  5.  Anamorphic  image. 


Fig.  6.  Fourier  synthesis  of  texture:  top  shows  texture  pattern; 
bottom,  texture  object. 


Fig.  7.  Plain  teapot  with  highlights. 
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( Xr , Yr,  Zr).  (Such  maps  will  be  shown  with  azimuthal 
angle  plotted  as  abscissa  and  polar  angle  plotted  as 
ordinate.) 

When  a reflection  direction  is  computed,  it  is  con 
verted  to  polar  coordinates  and  the  reflected  light  in 
tensity  for  that  direction  is  read  from  the  map.  This  is 
similar  to  the  technique  used  to  map  texture  onto  a 
surface,  except  that  the  reflection  direction,  instead  of 
parametric  surface  position,  determines  the  coordi 
nates  in  the  map.  Figure  8 shows  an  image  generated 
using  these  techniques. 

Use  of  the  surface  normal  alone  is  tantamount  to 
modeling  the  environment  on  an  infinitely  large  sphere. 
This  has  the  undesirable  effect  that  the  reflected  in 
tensity  at  all  silhouette  points  on  the  object  is  the  same, 
and  corresponds  to  the  intensity  of  the  environment 
model  diametrically  opposite  the  eye.  This  deficiency 
can  be  corrected  by  using  both  the  surface  normal 
and  surface  position  to  determine  what  part  of  the 
environment  map  is  reflected  in  a given  surface 
fragment. 

Combinations  of  Techniques 

The  techniques  described  above  for  simulating  tex- 
ture and  reflection  can  be  combined  to  produce  images 
of  objects  having  patterned  shiny  surfaces.  When  high 
lighting  is  combined  with  texture  mapping,  only  the 
component  from  the  real  light  source  should  be  scaled. 
This  models  the  highlight  as  being  specularly  reflected 
at  the  surface  and  not  being  affected  by  the  pigment 
within  it.  In  Figure  9 note  how  the  highlights  wash 
out  the  texture  pattern  underneath  them.  The  technique 
for  texture  mapping  actually  keys  the  texture  to  the 
surface  so  that  it  moves  with  the  object.  Some  other 
techniques,  which  essentially  apply  texture  to  the  2-D 
image,  do  not  have  this  property.  Note,  in  Figure  9, 
how  the  highlights  hardly  move  with  the  teapot,  whereas 
the  texture  does. 

Given  the  texture  mapping  technique  and  the  en- 
vironment reflecting  technique,  we  can  combine  them 


Fig.  9.  Textured  object  with  highlights,  two  orientations. 


Fig  10.  Highly  glazed  patterned  teapot. 
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to  produce  an  image  of  a highly  glazed  patterned  teapot, 
as  in  Figure  10. 


Resource  Requirements 

The  images  shown  in  this  paper  were  all  generated 
on  a PDP-11/45  computer  having  a 256K  byte  random 
access  frame  buffer  which  was  used  as  the  depth  buffer. 
The  main  routines  were  written  in  Fortran  and  the 
critical  parts  were  written  in  assembly  language.  The 
computation  time  of  the  extended  subdivision  algo- 
rithm is  roughly  proportional  to  the  area  covered  by 
visible  objects.  Images  of  nontextured  objects  of 
the  type  used  in  this  paper  take  about  25  minutes.  The 
addition  of  texture  or  reflection  increases  this  time  by 
about  10  percent.  All  images  have  a resolution  of 
512X512  picture  elements. 


Conclusions 


* 


By  refining  and  extending  Catmull’s  subdivision 
algorithm,  images  can  be  generated  having  a far  higher 
degree  of  naturalness  than  was  previously  possible. 

These  generalizations  result  in  improved  techniques 
for  generating  patterns  and  texture,  and  in  the  new 
capability  for  simulating  reflections. 
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AUDIO  PROCESSING 


At  the  end  of  September  1976  the  audio  processing  work 
being  pursued  under  t^  i s contract  became  separately  iunded 
under  ARPA  order  3301,  contract  U00 1 73-77-C-004 1 with  Naval 
Research  Labs.  The  small  amount  of  prepatory  work  done 
during  the  three  months  July  1976  - September  1976  is 
reported  as  part  of  the  first  semi-annual  technical  report 
under  that  contract:  Report  Number  UTEC-CSc-77-090 . 
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